#include "DSelector_cascade_mm.h"

void DSelector_cascade_mm::Init(TTree *locTree)
{
	// The Init() function is called when the selector needs to initialize a new tree or chain.
	// Typically here the branch addresses and branch pointers of the tree will be set.
	// Init() will be called many times when running on PROOF (once per file to be processed).

	//SET OUTPUT FILE NAME //can be overriden by user in PROOF
	dOutputFileName = "cascade_mm.root"; //"" for none
	dOutputTreeFileName = ""; //"" for none

	//DO THIS NEXT
	//Because this function gets called for each TTree in the TChain, we must be careful:
		//We need to re-initialize the tree interface & branch wrappers, but don't want to recreate histograms
	bool locInitializedPriorFlag = dInitializedFlag; //save whether have been initialized previously
	DSelector::Init(locTree); //This must be called to initialize wrappers for each new TTree
	//gDirectory now points to the output file with name dOutputFileName (if any)

	dBeamBunchPeriod = dAnalysisUtilities.Get_BeamBunchPeriod(Get_RunNumber());
	if(locInitializedPriorFlag)
		return; //have already created histograms, etc. below: exit

	//THEN THIS
	Get_ComboWrappers();
	dPreviousRunNumber = 0;

	/*********************************** EXAMPLE USER INITIALIZATION: ANALYSIS ACTIONS **********************************/

	//ANALYSIS ACTIONS: //Executed in order if added to dAnalysisActions
	//false/true below: use measured/kinfit data

	//PID
	dAnalysisActions.push_back(new DHistogramAction_ParticleID(dComboWrapper, false));
	dAnalysisActions.push_back(new DCutAction_dEdx(dComboWrapper, false, KPlus, SYS_CDC));
	dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 0.2, KPlus, SYS_BCAL));
	dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 0.15, KPlus, SYS_TOF));
	dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 0.2, KPlus, SYS_FCAL));
	//dAnalysisActions.push_back(new DCutAction_Kinematics(dComboWrapper, 0, KPlus, false, 0.0, 1.2));
	dAnalysisActions.push_back(new DHistogramAction_ParticleID(dComboWrapper, false, "PostPID"));

	//MASSES
	dAnalysisActions.push_back(new DHistogramAction_MissingMass(dComboWrapper, false, 3000, 0.5, 3.5));
	dAnalysisActions.push_back(new DHistogramAction_MissingMass(dComboWrapper, true, 3000, 0.5, 3.5, "KinFit"));
	dAnalysisActions.push_back(new DHistogramAction_MissingMass(dComboWrapper, true, 3000, 0.5, 3.5, "KinFit_Center"));
	(static_cast<DHistogramAction_MissingMass*>(dAnalysisActions.back()))->dBeamBunchRange = pair<int, int>(0, 0);
	dAnalysisActions.push_back(new DHistogramAction_MissingMass(dComboWrapper, true, 3000, 0.5, 3.5, "KinFit_Sideband"));
	(static_cast<DHistogramAction_MissingMass*>(dAnalysisActions.back()))->dBeamBunchRange = pair<int, int>(1, 1);

	//KINFIT RESULTS
	dAnalysisActions.push_back(new DHistogramAction_KinFitResults(dComboWrapper));
	dAnalysisActions.push_back(new DCutAction_KinFitFOM(dComboWrapper, 1.0E-4));

	//MASSES
	dAnalysisActions.push_back(new DHistogramAction_MissingMass(dComboWrapper, false, 3000, 0.5, 3.5, "PostKinFitCut"));
	dAnalysisActions.push_back(new DHistogramAction_MissingMass(dComboWrapper, true, 3000, 0.5, 3.5, "PostKinFitCut_KinFit"));
	dAnalysisActions.push_back(new DHistogramAction_MissingMass(dComboWrapper, true, 3000, 0.5, 3.5, "PostKinFitCut_KinFit_Center"));
	(static_cast<DHistogramAction_MissingMass*>(dAnalysisActions.back()))->dBeamBunchRange = pair<int, int>(0, 0);
	dAnalysisActions.push_back(new DHistogramAction_MissingMass(dComboWrapper, true, 3000, 0.5, 3.5, "PostKinFitCut_KinFit_Sideband"));
	(static_cast<DHistogramAction_MissingMass*>(dAnalysisActions.back()))->dBeamBunchRange = pair<int, int>(1, 1);

	//BEAM ENERGY
	dAnalysisActions.push_back(new DHistogramAction_BeamEnergy(dComboWrapper, false));
	//dAnalysisActions.push_back(new DCutAction_BeamEnergy(dComboWrapper, false, 8.4, 9.05));

	//KINEMATICS
	dAnalysisActions.push_back(new DHistogramAction_ParticleComboKinematics(dComboWrapper, false));

	//INITIALIZE ACTIONS
	//If you create any actions that you want to run manually (i.e. don't add to dAnalysisActions), be sure to initialize them here as well
	Initialize_Actions();

	/***************************************** CUSTOM HISTOGRAMS ****************************************/

	dMissMassOffKPluses_BestCombo_RFSignal = new TH2D("MissMassOffKPluses_BestCombo_RFSignal", ";Beam Energy (GeV);#gammap#rightarrowK^{#plus}K^{#plus} Missing Mass (GeV/c^{2})", 900, 3.0, 12.0, 1500, 0.5, 3.5);
	dMissMassOffKPluses_BestCombo_RFSideband = new TH2D("MissMassOffKPluses_BestCombo_RFSideband", ";Beam Energy (GeV);#gammap#rightarrowK^{#plus}K^{#plus} Missing Mass (GeV/c^{2})", 900, 3.0, 12.0, 1500, 0.5, 3.5);
	dMissMassOffKPluses_BestCombo_RFCut = new TH2D("MissMassOffKPluses_BestCombo_RFCut", ";Beam Energy (GeV);#gammap#rightarrowK^{#plus}K^{#plus} Missing Mass (GeV/c^{2})", 900, 3.0, 12.0, 1500, 0.5, 3.5);

	/***************************************** ADVANCED: CHOOSE BRANCHES TO READ ****************************************/

	//TO SAVE PROCESSING TIME
		//If you know you don't need all of the branches/data, but just a subset of it, you can speed things up
		//By default, for each event, the data is retrieved for all branches
		//If you know you only need data for some branches, you can skip grabbing data from the branches you don't need
		//Do this by doing something similar to the commented code below

	//dTreeInterface->Clear_GetEntryBranches(); //now get none
	//dTreeInterface->Register_GetEntryBranch("Proton__P4"); //manually set the branches you want
}

Bool_t DSelector_cascade_mm::Process(Long64_t locEntry)
{
	// The Process() function is called for each entry in the tree. The entry argument
	// specifies which entry in the currently loaded tree is to be processed.
	//
	// This function should contain the "body" of the analysis. It can contain
	// simple or elaborate selection criteria, run algorithms on the data
	// of the event and typically fill histograms.
	//
	// The processing can be stopped by calling Abort().
	// Use fStatus to set the return value of TTree::Process().
	// The return value is currently not used.

	//CALL THIS FIRST
	DSelector::Process(locEntry); //Gets the data from the tree for the entry
	//cout << "RUN " << Get_RunNumber() << ", EVENT " << Get_EventNumber() << endl;
	//TLorentzVector locProductionX4 = Get_X4_Production();

	/******************************************** GET POLARIZATION ORIENTATION ******************************************/

	//Only if the run number changes
	//RCDB environment must be setup in order for this to work! (Will return false otherwise)
	UInt_t locRunNumber = Get_RunNumber();
	if(locRunNumber != dPreviousRunNumber)
	{
		dIsPolarizedFlag = dAnalysisUtilities.Get_IsPolarizedBeam(locRunNumber, dIsPARAFlag);
		dPreviousRunNumber = locRunNumber;
	}

	/********************************************* SETUP UNIQUENESS TRACKING ********************************************/

	//ANALYSIS ACTIONS: Reset uniqueness tracking for each action
	//For any actions that you are executing manually, be sure to call Reset_NewEvent() on them here
	Reset_Actions_NewEvent();

	//PREVENT-DOUBLE COUNTING WHEN HISTOGRAMMING
		//Sometimes, some content is the exact same between one combo and the next
			//e.g. maybe two combos have different beam particles, but the same data for the final-state
		//When histogramming, you don't want to double-count when this happens: artificially inflates your signal (or background)
		//So, for each quantity you histogram, keep track of what particles you used (for a given combo)
		//Then for each combo, just compare to what you used before, and make sure it's unique

	//EXAMPLE 1: Particle-specific info:
	set<Int_t> locUsedSoFar_BeamEnergy; //Int_t: Unique ID for beam particles. set: easy to use, fast to search

	//EXAMPLE 2: Combo-specific info:
		//In general: Could have multiple particles with the same PID: Use a set of Int_t's
		//In general: Multiple PIDs, so multiple sets: Contain within a map
		//Multiple combos: Contain maps within a set (easier, faster to search)
	set<map<Particle_t, set<Int_t> > > locUsedSoFar_MissingMass;

	//INSERT USER ANALYSIS UNIQUENESS TRACKING HERE

	/************************************************* LOOP OVER COMBOS *************************************************/

	//Loop over combos
	double locBestBeamE = -1.0, locBestMissingMass_OffKPluses = -1.0;
	bool locRFSignalFlag = false;
	double locBestKinFitConfidenceLevel = -1.0;

	double locBestBeamE_RFCut = -1.0, locBestMissingMass_OffKPluses_RFCut = -1.0, locBestKinFitConfidenceLevel_RFCut = -1.0;
	for(UInt_t loc_i = 0; loc_i < Get_NumCombos(); ++loc_i)
	{
		//Set branch array indices for combo and all combo particles
		dComboWrapper->Set_ComboIndex(loc_i);

		// Is used to indicate when combos have been cut
		if(dComboWrapper->Get_IsComboCut()) // Is false when tree originally created
			continue; // Combo has been cut previously

		/********************************************** GET PARTICLE INDICES *********************************************/

		//Used for tracking uniqueness when filling histograms, and for determining unused particles

		//Step 0
		Int_t locBeamID = dComboBeamWrapper->Get_BeamID();
		Int_t locKPlus1TrackID = dKPlus1Wrapper->Get_TrackID();
		Int_t locKPlus2TrackID = dKPlus2Wrapper->Get_TrackID();

		/*********************************************** GET FOUR-MOMENTUM **********************************************/

		// Get P4's: //is kinfit if kinfit performed, else is measured
		//dTargetP4 is target p4
		//Step 0
		TLorentzVector locBeamP4 = dComboBeamWrapper->Get_P4();
		TLorentzVector locKPlus1P4 = dKPlus1Wrapper->Get_P4();
		TLorentzVector locKPlus2P4 = dKPlus2Wrapper->Get_P4();

		// Get Measured P4's:
		//Step 0
		TLorentzVector locBeamP4_Measured = dComboBeamWrapper->Get_P4_Measured();
		TLorentzVector locKPlus1P4_Measured = dKPlus1Wrapper->Get_P4_Measured();
		TLorentzVector locKPlus2P4_Measured = dKPlus2Wrapper->Get_P4_Measured();

		/********************************************* COMBINE FOUR-MOMENTUM ********************************************/

		// DO YOUR STUFF HERE

		// Combine 4-vectors
		TLorentzVector locMissingP4_Measured = locBeamP4_Measured + dTargetP4;
		locMissingP4_Measured -= locKPlus1P4_Measured + locKPlus2P4_Measured;

		TLorentzVector locMissingP4_OffKPluses = locBeamP4 + dTargetP4 - locKPlus1P4 - locKPlus2P4;
		double locMissingMass_OffKPluses = locMissingP4_OffKPluses.M();

		/******************************************** EXECUTE ANALYSIS ACTIONS *******************************************/

		// Loop through the analysis actions, executing them in order for the active particle combo
		if(!Execute_Actions()) //if the active combo fails a cut, IsComboCutFlag automatically set
			continue;

		//get beam/rf delta-t
		TLorentzVector locBeamX4 = dComboBeamWrapper->Get_X4();
		double locRFTime = dComboWrapper->Get_RFTime_Measured();
		double locBeamRFDeltaT = locBeamX4.T() - (locRFTime + (locBeamX4.Z() - dTargetCenter.Z())/29.9792458);

		double locKinFitConfidenceLevel = dComboWrapper->Get_ConfidenceLevel_KinFit();

		//RF CUT
		if(fabs(locBeamRFDeltaT) <= (0.5*dBeamBunchPeriod))
		{
			if(locKinFitConfidenceLevel <= locBestKinFitConfidenceLevel_RFCut)
				continue;

			locBestKinFitConfidenceLevel_RFCut = locKinFitConfidenceLevel;
			locBestBeamE_RFCut = locBeamP4.E();
			locBestMissingMass_OffKPluses_RFCut = locMissingMass_OffKPluses;
		}

		if(locKinFitConfidenceLevel <= locBestKinFitConfidenceLevel)
			continue;

		locBestKinFitConfidenceLevel = locKinFitConfidenceLevel;
		locBestBeamE = locBeamP4.E();
		locBestMissingMass_OffKPluses = locMissingMass_OffKPluses;

		//set flags
		locRFSignalFlag = fabs(locBeamRFDeltaT) <= (0.5*dBeamBunchPeriod);
	} // end of combo loop

	//FILL HISTOGRAMS: Num combos / events surviving actions
	Fill_NumCombosSurvivedHists();

	if(locBestBeamE > 0.0)
	{
		if(locRFSignalFlag)
			dMissMassOffKPluses_BestCombo_RFSignal->Fill(locBestBeamE, locBestMissingMass_OffKPluses);
		else
			dMissMassOffKPluses_BestCombo_RFSideband->Fill(locBestBeamE, locBestMissingMass_OffKPluses);
	}

	if(locBestBeamE_RFCut > 0.0)
		dMissMassOffKPluses_BestCombo_RFCut->Fill(locBestBeamE_RFCut, locBestMissingMass_OffKPluses_RFCut);

	/******************************************* LOOP OVER THROWN DATA (OPTIONAL) ***************************************/
/*
	//Thrown beam: just use directly
	if(dThrownBeam != NULL)
		double locEnergy = dThrownBeam->Get_P4().E();

	//Loop over throwns
	for(UInt_t loc_i = 0; loc_i < Get_NumThrown(); ++loc_i)
	{
		//Set branch array indices corresponding to this particle
		dThrownWrapper->Set_ArrayIndex(loc_i);

		//Do stuff with the wrapper here ...
	}
*/
	/****************************************** LOOP OVER OTHER ARRAYS (OPTIONAL) ***************************************/
/*
	//Loop over beam particles (note, only those appearing in combos are present)
	for(UInt_t loc_i = 0; loc_i < Get_NumBeam(); ++loc_i)
	{
		//Set branch array indices corresponding to this particle
		dBeamWrapper->Set_ArrayIndex(loc_i);

		//Do stuff with the wrapper here ...
	}

	//Loop over charged track hypotheses
	for(UInt_t loc_i = 0; loc_i < Get_NumChargedHypos(); ++loc_i)
	{
		//Set branch array indices corresponding to this particle
		dChargedHypoWrapper->Set_ArrayIndex(loc_i);

		//Do stuff with the wrapper here ...
	}

	//Loop over neutral particle hypotheses
	for(UInt_t loc_i = 0; loc_i < Get_NumNeutralHypos(); ++loc_i)
	{
		//Set branch array indices corresponding to this particle
		dNeutralHypoWrapper->Set_ArrayIndex(loc_i);

		//Do stuff with the wrapper here ...
	}
*/

	/************************************ EXAMPLE: FILL CLONE OF TTREE HERE WITH CUTS APPLIED ************************************/
/*
	Bool_t locIsEventCut = true;
	for(UInt_t loc_i = 0; loc_i < Get_NumCombos(); ++loc_i) {
		//Set branch array indices for combo and all combo particles
		dComboWrapper->Set_ComboIndex(loc_i);
		// Is used to indicate when combos have been cut
		if(dComboWrapper->Get_IsComboCut())
			continue;
		locIsEventCut = false; // At least one combo succeeded
		break;
	}
	if(!locIsEventCut && dOutputTreeFileName != "")
		FillOutputTree();
*/

	return kTRUE;
}

void DSelector_cascade_mm::Finalize(void)
{
	//Save anything to output here that you do not want to be in the default DSelector output ROOT file.

	//Otherwise, don't do anything else (especially if you are using PROOF).
		//If you are using PROOF, this function is called on each thread,
		//so anything you do will not have the combined information from the various threads.
		//Besides, it is best-practice to do post-processing (e.g. fitting) separately, in case there is a problem.

	//DO YOUR STUFF HERE

	//CALL THIS LAST
	DSelector::Finalize(); //Saves results to the output file
}