// $Id$
//
//    File: DCustomAction_p2pi_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_p2pi_hists.h"

void DCustomAction_p2pi_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();
		
		dM2pi = GetOrCreate_Histogram<TH1D>("M2pi", "#pi^{+} #pi^{-} Mass; Mass (GeV/c^{2})", 200, 0.0, 2.0);
		dMM_M2pi = GetOrCreate_Histogram<TH2D>("MM_M2pi", "MM off #pi^{+}#pi^{-} vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; MM", 200, 0.0, 2.0, 200, 0., 4.);
		dEgamma2pi = GetOrCreate_Histogram<TH1D>("Egamma2pi", "TAGGER photon energy; E_{gamma}", 320, 2., 10.);
		dMissingDeltaPT_DeltaPhi = GetOrCreate_Histogram<TH2D>("MissingDeltaPT_DeltaPhi", "Unused Track - MM: #Delta p_{T} vs #Delta #phi; #Delta #phi; #Delta p_{T}", 180, -180.0, 180.0, 100, -5., 5.);
                dMM2_DeltaPhi = GetOrCreate_Histogram<TH2D>("MM2_DeltaPhi", "MM^{2} vs #Delta #phi; #Delta #phi; MM^{2}", 180, -60.0, 60.0, 200, -0.5, 0.5);
		dMM_M2pi_Emeson = GetOrCreate_Histogram<TH2D>("MM_M2pi_Emeson", "MM off #pi^{+}#pi^{-} vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; MM", 200, 0.0, 2.0, 200, 0., 4.);
		
		// recoil and exclusive
		dMM_M2pi_Recoil = GetOrCreate_Histogram<TH2D>("MM_M2pi_Recoil", "Recoil MM off #pi^{+}#pi^{-} vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; MM", 200, 0.0, 2.0, 200, 0., 4.);
		dMM_M2pi_Exclusive = GetOrCreate_Histogram<TH2D>("MM_M2pi_Exclusive", "Exclusive MM off #pi^{+}#pi^{-} vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; MM", 200, 0.0, 2.0, 200, 0., 4.);
		dProton_dEdx_P_Exclusive = GetOrCreate_Histogram<TH2D>("Proton_dEdx_P_Exclusive","Exclusive: dE/dx vs p; p; dE/dx",200,0,2,500,0,5);
		
		dProton_P_Theta_Exclusive = GetOrCreate_Histogram<TH2D>("Proton_P_Theta_Exclusive","Exclusive: p vs #theta; #theta; p (GeV/c)",180,0,180,120,0,12);
		
		dDeltaE_M2pi_Exclusive = GetOrCreate_Histogram<TH2D>("dDeltaE_M2pi_Exclusive", "Exclusive #Delta E vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; #Delta E (tagger - tracks)", 200, 0.0, 2.0, 200, -5., 5.);
		dRecoilPT_M2pi = GetOrCreate_Histogram<TH2D>("dRecoilPT_M2pi", "Recoil P_{T} vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; Recoil P_{T}", 200, 0.0, 2.0, 200, 0., 2.);
		dRecoilPT_M2pi_Exclusive = GetOrCreate_Histogram<TH2D>("dRecoilPT_M2pi_Exclusive", "Exclusive Recoil P_{T} vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; Recoil P_{T}", 200, 0.0, 2.0, 200, 0., 2.);

		// MC histograms
		dMM_DeltaM2pi = GetOrCreate_Histogram<TH2D>("MM_DeltaM2pi", "MM off #pi^{+}#pi^{-} vs #Delta M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; MM", 200, -0.5, 0.5, 200, -4., 4.);
		dMM2_DeltaM2pi = GetOrCreate_Histogram<TH2D>("MM2_DeltaM2pi", "MM^{2} off #pi^{+}#pi^{-} vs #Delta M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; MM^{2}", 200, -0.5, 0.5, 200, -0.5, 0.5);
		dCosThetaCM_t = GetOrCreate_Histogram<TH2D>("CosThetaCM_t", "Mandelstam -t vs cos#theta_{CM}^{#rho}; cos#theta_{CM}^{#rho}; -t", 1000, -1., 1., 1000, 0., 5.);

	}
	japp->RootUnLock(); //RELEASE ROOT LOCK!!
}

bool DCustomAction_p2pi_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 DChargedTrack*> locUnusedChargedTracks;

	// should only have one reaction step
	const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(0);
	if(locParticleComboStep->Get_InitialParticleID() != Gamma)
		return true;
	
	// get beam photon time and select energy range
        const DKinematicData* locBeamPhoton = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle();
	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);

	// sum 2pi mass
	DLorentzVector locP4_2pi;
	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;
		locP4_2pi += locParticles[loc_i]->lorentzMomentum();
	}

	// loop over charged tracks to make list of unused
	for(size_t loc_j = 0; loc_j < locChargedTracks.size(); ++loc_j) {

		int nMatched = 0;

		// eliminate pi+/pi- from unused tracks
		for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) {
			
			// 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();
			if(locChargedTracks[loc_j]->Get_BestFOM()->candidateid == locChargedTrackHypothesis->candidateid)
				nMatched++;
		}
		
		if(nMatched==0)
			locUnusedChargedTracks.push_back(locChargedTracks[loc_j]);
	}

	// calculate missing mass
	DLorentzVector locMissingP4 = dAnalysisUtilities->Calc_MissingP4(locParticleCombo,Get_UseKinFitResultsFlag());

	// get 2pi mass from MC tracks
	vector<const DMCThrown*> locMCThrown;
	locEventLoop->Get(locMCThrown);
	DLorentzVector locMCP4_2pi;
	if(!locMCThrown.empty()){
		if(!Get_UseKinFitResultsFlag()) {
			//cout<<"Found good event "<<locEventLoop->GetJEvent().GetEventNumber()<<": M_2pi="<<locP4_2pi.M()<<" MM="<<locMissingP4.M()<<endl;
			
			int nPiP = 0;
			int nPiM = 0;
			for(uint i=0; i<locMCThrown.size(); i++){
				if(locMCThrown[i]->type == 8){
					locMCP4_2pi += locMCThrown[i]->lorentzMomentum();
					nPiP++;
				}
				if(locMCThrown[i]->type == 9){
					locMCP4_2pi += locMCThrown[i]->lorentzMomentum();
					nPiM++;
				}
			}
		}
	}

	japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
	{
		// Fill any histograms here
		dM2pi->Fill(locP4_2pi.M());
		dMM_M2pi->Fill(locP4_2pi.M(), locMissingP4.M());
		dEgamma2pi->Fill(locBeamPhotonEnergy);
		if(!locMCThrown.empty()) 
			dMM_DeltaM2pi->Fill((locP4_2pi.M()-locMCP4_2pi.M()), locMissingP4.M());
		
		double missingMomMag = 0.3;

		if(locMissingP4.Vect().Perp() > missingMomMag) 
			dMM_M2pi_Recoil->Fill(locP4_2pi.M(), locMissingP4.M());
		dRecoilPT_M2pi->Fill(locP4_2pi.M(),locMissingP4.Vect().Perp());

		// correlate best hypotheses with 2pi system
		for(size_t loc_j = 0; loc_j < locUnusedChargedTracks.size(); ++loc_j) {
			const DChargedTrackHypothesis* locChargedTrackHypothesis = locUnusedChargedTracks[loc_j]->Get_Hypothesis(Proton);
			DLorentzVector locUnusedP4;
			
			// if no proton hypothesis, need to make ad-hoc version with correct mass
			if(locChargedTrackHypothesis != NULL) 
				locUnusedP4 = locChargedTrackHypothesis->lorentzMomentum();
			else { 
				locChargedTrackHypothesis = locUnusedChargedTracks[loc_j]->Get_BestTrackingFOM();
				DLorentzVector locUnusedP4 = locChargedTrackHypothesis->lorentzMomentum();
				locUnusedP4.SetVectM(locUnusedP4.Vect(), 0.938);
			}
	
                        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
			DLorentzVector locMissingTotalP4(0,0,locBeamPhotonEnergy,locBeamPhotonEnergy+0.938);
			locMissingTotalP4 -= locUnusedP4;
			locMissingTotalP4 -= locP4_2pi;
			dMM2_DeltaPhi->Fill(deltaPhi, locMissingTotalP4.M2());
			if(!locMCThrown.empty()) 
				dMM2_DeltaM2pi->Fill((locP4_2pi.M()-locMCP4_2pi.M()), locMissingTotalP4.M2());
			
			// plot exclusive results
			if(fabs(deltaPhi) < 5 && fabs(locMissingTotalP4.M2()) < 0.1){
				dMM_M2pi_Exclusive->Fill(locP4_2pi.M(), locMissingP4.M());
				dRecoilPT_M2pi_Exclusive->Fill(locP4_2pi.M(),locMissingP4.Vect().Perp());

				if(locMissingP4.Vect().Perp() > missingMomMag) {	
					double dEdx = locChargedTrackHypothesis->dEdx()*1e3;
					if(!locMCThrown.empty()) dEdx *= 1e2; // scale MC dE/dx
					dProton_dEdx_P_Exclusive->Fill(locUnusedP4.Vect().Mag(), dEdx);
					dProton_P_Theta_Exclusive->Fill(locUnusedP4.Vect().Theta()*180/TMath::Pi(), locUnusedP4.Vect().Mag());
					
					dDeltaE_M2pi_Exclusive->Fill(locP4_2pi.M(),locMissingTotalP4.E());
				}
			}
                }

		// 2pi energy cut (Dave M. suggestion)
		if(locP4_2pi.E() > 5) {
			dMM_M2pi_Emeson->Fill(locP4_2pi.M(), locMissingP4.M());
		}

		// 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 loc2piP4CMFrame = locP4_2pi;
                locGammaP4CMFrame.Boost(locBoostVector);
                loc2piP4CMFrame.Boost(locBoostVector);

                double locThetaCM = locGammaP4CMFrame.Vect().Angle(loc2piP4CMFrame.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!)
}