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

void DCustomAction_sdobbs_pi0_dalitz::Initialize(JEventLoop* locEventLoop)
{

    /*  -- disable for analysis launch

	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();

		dEgamma = GetOrCreate_Histogram<TH1I>("Egamma", "TAGGER photon energy; E_{#gamma}", 400, 0., 12.);
	
		dmass_ee_v_masseeg = GetOrCreate_Histogram<TH2I>("mass e^{+}e^{-}","mass e^{+}e^{-}",400,0.,2.,400,0.,2.);
		ebyplogic_mee_v_meeg = GetOrCreate_Histogram<TH2I>("mass e^{+}e^{-} logic1","mass e^{+}e^{-} logic1",400,0.,2.,400,0.,2.);
		dmassPi0_ee_cut = GetOrCreate_Histogram<TH1I>("mass e^{+}e^{-}#gamma post_mass_e^{+}e^{-} cut","mass e^{+}e^{-}#gamma post_mass_e^{+}e^{-} cut",400,0.,2.);
		ebyplogic_both_mee_v_meeg = GetOrCreate_Histogram<TH2I>("mass e^{+}e^{-} logic2","mass e^{+}e^{-} logic2",400,0.,2.,400,0.,2.);
		meeg_v_pmag = GetOrCreate_Histogram<TH2I>("meeg v pmag","meeg v pmag",200,0,5,200,0,2);
		meeg_v_shwrE = GetOrCreate_Histogram<TH2I>("meeg v Eshwr","meeg v Eshwr",200,0,2,200,0,2);
		BCAL_Ee_minus_pe_vs_pe = GetOrCreate_Histogram<TH2I>("BCAL electron Ee-pe v pe","BCAL electron Ee-pe v pe",200,0,2,200,-2,2);
		BCAL_Ee_plus_pe_vs_pe = GetOrCreate_Histogram<TH2I>("BCAL positron Ee-pe v pe","BCAL positron Ee-pe v pe",200,0,2,400,-2,2);
		FCAL_Ee_minus_pe_vs_pe = GetOrCreate_Histogram<TH2I>("FCAL electron Ee-pe v pe","FCAL electron Ee-pe v pe",200,0,2,200,-2,2);
		FCAL_Ee_plus_pe_vs_pe = GetOrCreate_Histogram<TH2I>("FCAL positron Ee-pe v pe","FCAL positron Ee-pe v pe",200,0,2,200,-2,2);

		ebyp = GetOrCreate_Histogram<TH1I>("ebyp","ebyp",200,0,2);
		ebyp_BCAL_plus = GetOrCreate_Histogram<TH1I>("ebyp BCAL plus","ebyp BCAL plus",200,0,2);
		ebyp_BCAL_neg = GetOrCreate_Histogram<TH1I>("ebyp BCAL neg","ebyp BCAL neg",200,0,2);
		ebyp_FCAL_plus = GetOrCreate_Histogram<TH1I>("ebyp FCAL plus","ebyp FCAL plus",200,0,2);
		ebyp_FCAL_neg = GetOrCreate_Histogram<TH1I>("ebyp FCAL neg","ebyp FCAL neg",200,0,2);


		dMM2_MPi0 = GetOrCreate_Histogram<TH2I>("MM2_MPi0", "MM^{2} off e^{+}e^{-}#gamma vs M_{e^{+}e^{-}#gamma}; M_{e^{+}e^{-}#gamma}; MM^{2}", 200, 0.0, 2.0, 200, -1., 1.);
		dProton_dEdx_P = GetOrCreate_Histogram<TH2I>("Proton_dEdx_P","dE/dx vs p; p; dE/dx",200,0,2,500,0,5);
		dProton_P_Theta = GetOrCreate_Histogram<TH2I>("Proton_P_Theta","p vs #theta; #theta; p (GeV/c)",180,0,180,120,0,12);
		
	}
	japp->RootUnLock(); //RELEASE ROOT LOCK!!

    */
}

bool DCustomAction_sdobbs_pi0_dalitz::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
{
    return true; // disable for analysis launch

	const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(0);

	const DDetectorMatches* locDetectorMatches = NULL;
	locEventLoop->GetSingle(locDetectorMatches);

	vector<const DBCALShower*> locBCALShowers;
	locEventLoop->Get(locBCALShowers);

	vector<const DFCALShower*> locFCALShowers;
	locEventLoop->Get(locFCALShowers);
		
	vector<const DFCALCluster*> locFCALClusters;
	locEventLoop->Get(locFCALClusters);

	// get beam photon energy and final state particles
        const DKinematicData* locBeamPhoton = NULL;
        deque<const DKinematicData*> locParticles;
        if(!Get_UseKinFitResultsFlag()) { //measured
		locBeamPhoton = locParticleComboStep->Get_InitialParticle_Measured();
                locParticleComboStep->Get_FinalParticles_Measured(locParticles);
	}
	else {
		locBeamPhoton = locParticleComboStep->Get_InitialParticle();
		locParticleComboStep->Get_FinalParticles(locParticles);
	}
        double locBeamPhotonEnergy = locBeamPhoton->energy();

	DLorentzVector locSumInitP4;
	locSumInitP4.SetXYZM(0, 0, 0, 0.938);
    locSumInitP4 += locBeamPhoton->lorentzMomentum();

	// calculate missing mass
	DLorentzVector locMissingP4 = dAnalysisUtilities->Calc_MissingP4(locParticleCombo, Get_UseKinFitResultsFlag());
	
	// calculate 3pi mass
	DLorentzVector locPi0P4 = dAnalysisUtilities->Calc_FinalStateP4(locParticleCombo, 1, Get_UseKinFitResultsFlag());
	DLorentzVector locElectronP4 = locParticleCombo->Get_ParticleComboStep(1)->Get_FinalParticle_Measured(0)->lorentzMomentum();
	DLorentzVector locPositronP4 = locParticleCombo->Get_ParticleComboStep(1)->Get_FinalParticle_Measured(1)->lorentzMomentum();
	double locgammaE = locParticleCombo->Get_ParticleComboStep(1)->Get_FinalParticle_Measured(2)->energy();

	// reconstructed proton variables
	double dEdx = 0.;
	DLorentzVector locProtonP4;
	const DChargedTrack* locChargedTrack = static_cast<const DChargedTrack*>(locParticleComboStep->Get_FinalParticle_SourceObject(1));
	const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTrack->Get_Hypothesis(Proton);
	dEdx = locChargedTrackHypothesis->dEdx()*1e6;
	locProtonP4 = locChargedTrackHypothesis->lorentzMomentum();	

	double dEdxCut = 2.2;

	DVector3 electron_mom = locParticleCombo->Get_ParticleComboStep(1)->Get_FinalParticle_Measured(0)->momentum();
	DVector3 positron_mom = locParticleCombo->Get_ParticleComboStep(1)->Get_FinalParticle_Measured(1)->momentum();	

	double BCAL_ebyp_e_min = 0.;
	double BCAL_ebyp_e_pos = 0.;
	double FCAL_ebyp_e_min = 0.;
	double FCAL_ebyp_e_pos = 0.;
	int ebyp_logic = 0;
	int ebyp_eminus_logic = 0;
	int ebyp_eplus_logic = 0;
	
	for(unsigned int i = 0 ; i < locBCALShowers.size(); i++){
		if (locDetectorMatches->Get_IsMatchedToTrack(locBCALShowers[i])){
			BCAL_ebyp_e_min = locBCALShowers[i]->E_raw/electron_mom.Mag();
			ebyp->Fill(BCAL_ebyp_e_min);
			ebyp_BCAL_neg->Fill(BCAL_ebyp_e_min);
			if(locPi0P4.M()>0.1&&locPi0P4.M()<0.16&&(locElectronP4+locPositronP4).M()<.05) BCAL_Ee_minus_pe_vs_pe->Fill(electron_mom.Mag(),(locBCALShowers[i]->E_raw-electron_mom.Mag()));
			BCAL_ebyp_e_pos = locBCALShowers[i]->E_raw/positron_mom.Mag();
			ebyp->Fill(BCAL_ebyp_e_pos);
			ebyp_BCAL_plus->Fill(BCAL_ebyp_e_pos);
			if(locPi0P4.M()>0.1&&locPi0P4.M()<0.16&&(locElectronP4+locPositronP4).M()<.05) BCAL_Ee_plus_pe_vs_pe->Fill(positron_mom.Mag(),(locBCALShowers[i]->E_raw-electron_mom.Mag()));
			if( (BCAL_ebyp_e_min > 0.8 && BCAL_ebyp_e_min < 1.2) || (BCAL_ebyp_e_pos>0.8 && BCAL_ebyp_e_pos<1.2) ) ebyp_logic=1;
			if(BCAL_ebyp_e_min > 0.8 && BCAL_ebyp_e_min < 1.2) ebyp_eminus_logic =1;
			if(BCAL_ebyp_e_pos>0.8 && BCAL_ebyp_e_pos<1.2) ebyp_eplus_logic=1;
		}
	}
			  
	for(unsigned int i = 0 ; i < locFCALShowers.size(); i++){
		if (locDetectorMatches->Get_IsMatchedToTrack(locFCALShowers[i])){
			const DFCALShower *s1 = locFCALShowers[i];
			vector<const DFCALCluster*> associated_clusters;
			s1->Get(associated_clusters);
			FCAL_ebyp_e_min = s1->getEnergy()/electron_mom.Mag();
			ebyp->Fill(FCAL_ebyp_e_min);
			ebyp_FCAL_neg->Fill(FCAL_ebyp_e_min);
			if(locPi0P4.M()>0.1&&locPi0P4.M()<0.16&&(locElectronP4+locPositronP4).M()<.05) FCAL_Ee_minus_pe_vs_pe->Fill(electron_mom.Mag(),(s1->getEnergy()-electron_mom.Mag()));
			FCAL_ebyp_e_pos = s1->getEnergy()/positron_mom.Mag();
			ebyp->Fill(FCAL_ebyp_e_pos);
			ebyp_FCAL_plus->Fill(FCAL_ebyp_e_pos);
			 if(locPi0P4.M()>0.1&&locPi0P4.M()<0.16&&(locElectronP4+locPositronP4).M()<.05) FCAL_Ee_plus_pe_vs_pe->Fill(positron_mom.Mag(),(s1->getEnergy()-electron_mom.Mag()));
			if( (FCAL_ebyp_e_min > 0.8 && FCAL_ebyp_e_min < 1.2) || (FCAL_ebyp_e_pos>0.8 && FCAL_ebyp_e_pos<1.2) ) ebyp_logic=1;
			if(FCAL_ebyp_e_min > 0.8 && FCAL_ebyp_e_min < 1.2) ebyp_eminus_logic=1;
			if(FCAL_ebyp_e_pos>0.8 && FCAL_ebyp_e_pos<1.2) ebyp_eplus_logic=1;
		}
	}

	japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
	{
		// Fill histograms here
		dEgamma->Fill(locBeamPhotonEnergy);

		dmass_ee_v_masseeg->Fill((locElectronP4+locPositronP4).M(),locPi0P4.M());

		if((locElectronP4+locPositronP4).M()<0.075) dmassPi0_ee_cut->Fill(locPi0P4.M());

		dMM2_MPi0->Fill(locPi0P4.M(), locMissingP4.M2());

		if(ebyp_logic==1) ebyplogic_mee_v_meeg->Fill((locElectronP4+locPositronP4).M(),locPi0P4.M());
		if(ebyp_eminus_logic==1 && ebyp_eplus_logic==1) ebyplogic_both_mee_v_meeg->Fill((locElectronP4+locPositronP4).M(),locPi0P4.M());
		if(ebyp_eminus_logic==1&&ebyp_eplus_logic==1&&(locElectronP4+locPositronP4).M()<0.13) meeg_v_pmag->Fill(electron_mom.Mag(),locPi0P4.M()); 
		if(ebyp_eminus_logic==1&&ebyp_eplus_logic==1&&(locElectronP4+locPositronP4).M()<0.13) meeg_v_shwrE->Fill(locgammaE,locPi0P4.M());

		double locDeltaPhi = (locProtonP4.Phi() - locPi0P4.Phi())*180./TMath::Pi();
		if(locDeltaPhi > 360.) locDeltaPhi -= 360.;
		if(locDeltaPhi < 0.) locDeltaPhi += 360.;
		
		if(fabs(locMissingP4.M2()) < 0.05 && fabs(locMissingP4.E()) < 0.5) {
			
			dProton_dEdx_P->Fill(locProtonP4.Vect().Mag(), dEdx);
			dProton_P_Theta->Fill(locProtonP4.Vect().Theta()*180/TMath::Pi(), locProtonP4.Vect().Mag());
			
		}
	}
	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!)
}