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

//static const double kK0Mass = 0.497614;
//static const double kPhiMass = 1.020;

void DCustomAction_jpsi_hists::Initialize(JEventLoop* locEventLoop)
{

	// check if a particle is missing
	auto locMissingPIDs = Get_Reaction()->Get_MissingPIDs();
	dMissingPID = locMissingPIDs.empty() ? Unknown : locMissingPIDs[0];

	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.);
        dMee = GetOrCreate_Histogram<TH1I>("Mee", "; M(e^{+}e^{-}), GeV", 200, 2., 4.);
        dMeep = GetOrCreate_Histogram<TH1I>("Meep", "; M(e^{+}e^{-}p), GeV", 200, 3., 5.);

        dBCAL_EP_pos = GetOrCreate_Histogram<TH1D>("BCAL_EP_pos", "BCAL E/p (q+)", 100, 0., 2.);   
        dBCAL_EP_p_pos = GetOrCreate_Histogram<TH2D>("BCAL_EP_p_pos", "BCAL p v. E/p (q+)", 100, 0., 8., 100, 0., 2.); 
        dBCAL_EP_neg = GetOrCreate_Histogram<TH1D>("BCAL_EP_neg", "BCAL E/p (q-)", 100, 0., 2.);   
        dBCAL_EP_p_neg = GetOrCreate_Histogram<TH2D>("BCAL_EP_p_neg", "BCAL p v. E/p (q-)", 100, 0., 8., 100, 0., 2.); 
        dBCAL_EP_proton = GetOrCreate_Histogram<TH1D>("BCAL_EP_proton", "BCAL E/p (proton)", 100, 0., 2.);   
        dBCAL_EP_p_proton = GetOrCreate_Histogram<TH2D>("BCAL_EP_p_proton", "BCAL p v. E/p (proton)", 100, 0., 8., 100, 0., 2.); 

        dFCAL_EP_pos = GetOrCreate_Histogram<TH1D>("FCAL_EP_pos", "FCAL E/p (q+)", 100, 0., 2.);   
        dFCAL_EP_p_pos = GetOrCreate_Histogram<TH2D>("FCAL_EP_p_pos", "FCAL p v. E/p (q+)", 100, 0., 8., 100, 0., 2.); 
        dFCAL_EP_neg = GetOrCreate_Histogram<TH1D>("FCAL_EP_neg", "FCAL E/p (q-)", 100, 0., 2.);   
        dFCAL_EP_p_neg = GetOrCreate_Histogram<TH2D>("FCAL_EP_p_neg", "FCAL p v. E/p (q-)", 100, 0., 8., 100, 0., 2.); 
        dFCAL_EP_proton = GetOrCreate_Histogram<TH1D>("FCAL_EP_proton", "FCAL E/p (proton)", 100, 0., 2.);   
        dFCAL_EP_p_proton = GetOrCreate_Histogram<TH2D>("FCAL_EP_p_proton", "FCAL p v. E/p (proton)", 100, 0., 8., 100, 0., 2.); 
	}
	japp->RootUnLock(); //RELEASE ROOT LOCK!!

}

bool DCustomAction_jpsi_hists::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
{
  const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(1);

  // get beam photon energy and final state particles
  const DKinematicData* locBeamPhoton = NULL;
  vector<const DKinematicData*> locParticles;
  if(!Get_UseKinFitResultsFlag()) { //measured
    locBeamPhoton = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle_Measured();
    locParticles = locParticleComboStep->Get_FinalParticles_Measured();
  }
  else {
    locBeamPhoton = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle();
    locParticles = locParticleComboStep->Get_FinalParticles();
  }
  double locBeamPhotonEnergy = locBeamPhoton->energy();
	
	// cut on tagger energy REMOVED
	//if( (locBeamPhotonEnergy < dEbeamMin) || (locBeamPhotonEnergy > dEbeamMax) ) 
	//	return true;

    // make E/p plots
    vector<const DParticleID *> dParticleID_algos;
    locEventLoop->Get(dParticleID_algos);
    if(dParticleID_algos.size()<1){
        _DBG_<<"Unable to get a DParticleID object! NO PID will be done!"<<endl;
        return RESOURCE_UNAVAILABLE;
    }
    const DParticleID *dParticleID = dParticleID_algos[0];

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

    // -- electron 
    {
        const DChargedTrack* locChargedTrack = static_cast<const DChargedTrack*>(locParticleComboStep->Get_FinalParticle_SourceObject(0));
        if(locChargedTrack != NULL) {
        const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTrack->Get_Hypothesis(PiMinus);
        if(locChargedTrackHypothesis != NULL) {
        const DTrackTimeBased* locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased();
        if(locTrackTimeBased != NULL) {
        double locP = locTrackTimeBased->momentum().Mag();

        shared_ptr<const DFCALShowerMatchParams> locFCALShowerMatchParams;
        shared_ptr<const DBCALShowerMatchParams> locBCALShowerMatchParams;
        bool foundFCAL = dParticleID->Get_BestFCALMatchParams(locTrackTimeBased, locDetectorMatches, locFCALShowerMatchParams);
        bool foundBCAL = dParticleID->Get_BestBCALMatchParams(locTrackTimeBased, locDetectorMatches, locBCALShowerMatchParams);

        if(foundFCAL){
            const DFCALShower* locFCALShower = dynamic_cast<const DFCALShower*>(locFCALShowerMatchParams->dFCALShower); //dShowerObject
            double locFCALEnergy = locFCALShower->getEnergy();

            Lock_Action(); //ACQUIRE ROOT LOCK!!
            {
                double FCAL_EP = locFCALEnergy/locP;
                dFCAL_EP_neg->Fill(FCAL_EP);
                dFCAL_EP_p_neg->Fill(locP,FCAL_EP);
            }
            Unlock_Action(); //RELEASE ROOT LOCK!!
        }

        if(foundBCAL){
            const DBCALShower* locBCALShower = dynamic_cast<const DBCALShower*>(locBCALShowerMatchParams->dBCALShower); //dShowerObject
            double locBCALEnergy = locBCALShower->E;
            
            Lock_Action(); //ACQUIRE ROOT LOCK!!
            {
                double BCAL_EP = locBCALEnergy/locP;                            
                dBCAL_EP_neg->Fill(BCAL_EP);
                dBCAL_EP_p_neg->Fill(locP,BCAL_EP);
            }
            Unlock_Action(); //RELEASE ROOT LOCK!!
        }
        }
        }
        }
    }

    // -- positron
    {
        const DChargedTrack* locChargedTrack = static_cast<const DChargedTrack*>(locParticleComboStep->Get_FinalParticle_SourceObject(1));
        if(locChargedTrack != NULL) {
        const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTrack->Get_Hypothesis(PiPlus);
        if(locChargedTrackHypothesis != NULL) {
        const DTrackTimeBased* locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased();
        if(locTrackTimeBased != NULL) {
        double locP = locTrackTimeBased->momentum().Mag();

        shared_ptr<const DFCALShowerMatchParams> locFCALShowerMatchParams;
        shared_ptr<const DBCALShowerMatchParams> locBCALShowerMatchParams;
        bool foundFCAL = dParticleID->Get_BestFCALMatchParams(locTrackTimeBased, locDetectorMatches, locFCALShowerMatchParams);
        bool foundBCAL = dParticleID->Get_BestBCALMatchParams(locTrackTimeBased, locDetectorMatches, locBCALShowerMatchParams);

        if(foundFCAL){
            const DFCALShower* locFCALShower = dynamic_cast<const DFCALShower*>(locFCALShowerMatchParams->dFCALShower); //dShowerObject
            double locFCALEnergy = locFCALShower->getEnergy();

            Lock_Action(); //ACQUIRE ROOT LOCK!!
            {
                double FCAL_EP = locFCALEnergy/locP;
                dFCAL_EP_pos->Fill(FCAL_EP);
                dFCAL_EP_p_pos->Fill(locP,FCAL_EP);
            }
            Unlock_Action(); //RELEASE ROOT LOCK!!
        }

        if(foundBCAL){
            const DBCALShower* locBCALShower = dynamic_cast<const DBCALShower*>(locBCALShowerMatchParams->dBCALShower); //dShowerObject
            double locBCALEnergy = locBCALShower->E;
            
            Lock_Action(); //ACQUIRE ROOT LOCK!!
            {
                double BCAL_EP = locBCALEnergy/locP;                            
                dBCAL_EP_pos->Fill(BCAL_EP);
                dBCAL_EP_p_pos->Fill(locP,BCAL_EP);
            }
            Unlock_Action(); //RELEASE ROOT LOCK!!
        }
        }
        }
        }
    }


    // -- proton
    if(dMissingPID != Proton) {
        const DChargedTrack* locChargedTrack = static_cast<const DChargedTrack*>(locParticleCombo->Get_ParticleComboStep(0)->Get_FinalParticle_SourceObject(1));
        if(locChargedTrack != NULL) {
        const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTrack->Get_Hypothesis(Proton);
        if(locChargedTrackHypothesis != NULL) {
        const DTrackTimeBased* locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased();
        if(locTrackTimeBased != NULL) {
        double locP = locTrackTimeBased->momentum().Mag();

        shared_ptr<const DFCALShowerMatchParams> locFCALShowerMatchParams;
        shared_ptr<const DBCALShowerMatchParams> locBCALShowerMatchParams;
        bool foundFCAL = dParticleID->Get_BestFCALMatchParams(locTrackTimeBased, locDetectorMatches, locFCALShowerMatchParams);
        bool foundBCAL = dParticleID->Get_BestBCALMatchParams(locTrackTimeBased, locDetectorMatches, locBCALShowerMatchParams);

        if(foundFCAL){
            const DFCALShower* locFCALShower = dynamic_cast<const DFCALShower*>(locFCALShowerMatchParams->dFCALShower); //dShowerObject
            double locFCALEnergy = locFCALShower->getEnergy();

            Lock_Action(); //ACQUIRE ROOT LOCK!!
            {
                double FCAL_EP = locFCALEnergy/locP;
                dFCAL_EP_proton->Fill(FCAL_EP);
                dFCAL_EP_p_proton->Fill(locP,FCAL_EP);
            }
            Unlock_Action(); //RELEASE ROOT LOCK!!
        }

        if(foundBCAL){
            const DBCALShower* locBCALShower = dynamic_cast<const DBCALShower*>(locBCALShowerMatchParams->dBCALShower); //dShowerObject
            double locBCALEnergy = locBCALShower->E;
            
            Lock_Action(); //ACQUIRE ROOT LOCK!!
            {
                double BCAL_EP = locBCALEnergy/locP;                            
                dBCAL_EP_proton->Fill(BCAL_EP);
                dBCAL_EP_p_proton->Fill(locP,BCAL_EP);
            }
            Unlock_Action(); //RELEASE ROOT LOCK!!
        }
        }
        }
        }
    }


    // calc some 4 vectors
    //DLorentzVector loc_electron = locParticles[0]->lorentzMomentum();BUGGED ATM


	Lock_Action();
	{
	  dEgamma->Fill(locBeamPhotonEnergy);
	  //dMee->Fill(loc_electron.M()); BUGGED ATM
	}
	Unlock_Action(); //ACQUIRE ROOT LOCK!!


    /*
	DLorentzVector locP4_ee = locParticles[0]->lorentzMomentum() + locParticles[1]->lorentzMomentum();
    DLorentzVector locP4_eep;
    if(dMissingPID != Proton) {   // if we've reconstructed the proton
        locP4_eep = locP4_ee + locParticles[2]->lorentzMomentum();
    } else {  // if we're looking at a missing proton situation
        // calculate missing P4
        DLorentzVector locMissingP4 = dAnalysisUtilities->Calc_MissingP4(locParticleCombo,Get_UseKinFitResultsFlag());

        locP4_eep = locP4_ee + locMissingP4;
    }

	Lock_Action(); //ACQUIRE ROOT LOCK!!
	{
		// Fill histograms here
        dMee->Fill( locP4_ee.M() );

        if( (locP4_ee.M() > 2.7) && (locP4_ee.M() < 3.3) ) 
            dMeep->Fill( locP4_eep.M() );
    }
	Unlock_Action(); //ACQUIRE ROOT LOCK!!
    */



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

	// reconstructed proton variables
	double dEdx = 0.;
	double CDC_dEdx = -1.;
	double FDC_dEdx = -1.;
	double CDC_dEdx_pi1 = -1.;
	double FDC_dEdx_pi1 = -1.;
	double CDC_dEdx_pi2 = -1.;
	double FDC_dEdx_pi2 = -1.;
	//DLorentzVector locProtonP4;
	const DChargedTrack* locChargedTrack = static_cast<const DChargedTrack*>(locParticleComboStep->Get_FinalParticle_SourceObject(2));
	const DChargedTrack* locChargedTrack_pi1 = static_cast<const DChargedTrack*>(locParticleComboStep->Get_FinalParticle_SourceObject(0));
	const DChargedTrack* locChargedTrack_pi2 = static_cast<const DChargedTrack*>(locParticleComboStep->Get_FinalParticle_SourceObject(1));
	const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTrack->Get_Hypothesis(Proton);
	const DChargedTrackHypothesis* locChargedTrackHypothesis_pi1 = locChargedTrack_pi1->Get_Hypothesis(PiPlus);
	const DChargedTrackHypothesis* locChargedTrackHypothesis_pi2 = locChargedTrack_pi2->Get_Hypothesis(PiMinus);
	dEdx = locChargedTrackHypothesis->dEdx()*1e6;
	//dEdx = locChargedTrackHypothesis->dEdx()*1.e3;
	DLorentzVector locProtonP4 = locChargedTrackHypothesis->lorentzMomentum();	
	DLorentzVector locPion1P4 = locChargedTrackHypothesis_pi1->lorentzMomentum();	
	DLorentzVector locPion1P2 = locChargedTrackHypothesis_pi2->lorentzMomentum();	
     const DTrackTimeBased* locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased();
	const DTrackTimeBased* locTrackTimeBased_pi1 = locChargedTrackHypothesis_pi1->Get_TrackTimeBased();
	const DTrackTimeBased* locTrackTimeBased_pi2 = locChargedTrackHypothesis_pi2->Get_TrackTimeBased();
	if(locTrackTimeBased != NULL) {
		CDC_dEdx = locTrackTimeBased->ddEdx_CDC*1e6;
		FDC_dEdx = locTrackTimeBased->ddEdx_FDC*1e6;
	}
	if(locTrackTimeBased_pi1 != NULL) {
		CDC_dEdx_pi1 = locTrackTimeBased->ddEdx_CDC*1e6;
		FDC_dEdx_pi1 = locTrackTimeBased->ddEdx_FDC*1e6;
	}
	if(locTrackTimeBased_pi2 != NULL) {
		CDC_dEdx_pi2 = locTrackTimeBased->ddEdx_CDC*1e6;
		FDC_dEdx_pi2 = locTrackTimeBased->ddEdx_FDC*1e6;
	}

	Lock_Action(); //ACQUIRE ROOT LOCK!!
	{
		// Fill histograms here
		dEgamma->Fill(locBeamPhotonEnergy);
		dM2pi->Fill(locP4_2pi.M());
		dMM->Fill(locMissingP4.M());
		if(fabs(locP4_2pi.M() - kK0Mass)  < 0.02)
			dMM_M2picut->Fill(locMissingP4.M());

		dMM_M2pi->Fill(locP4_2pi.M(), locMissingP4.M());
		dDeltaE_M2pi->Fill(locP4_2pi.M(),locMissingP4.E());
			
		dVertexDeltaZ_M2pi->Fill(locP4_2pi.M(), locDiplacedVertex.Z());
		dVertexRho_M2pi->Fill(locP4_2pi.M(), locDiplacedVertex.Perp());
		dVertexDeltaMag_M2pi->Fill(locP4_2pi.M(), locDiplacedVertex.Mag());
		
        DLorentzVector locPhiP4 = locP4_2pi + locMissingP4;

        // ---
		if( (fabs(locMissingP4.M2() - kK0Mass*kK0Mass) < 0.05) 
		    && (fabs(locP4_2pi.M() - kK0Mass)  < 0.1) ) {
            dKL_Pmag->Fill( locMissingP4.Vect().Mag() );

			dProton_dEdx_P->Fill(locProtonP4.Vect().Mag(), dEdx);
			dProton_CDCdEdx_P->Fill(locProtonP4.Vect().Mag(), CDC_dEdx);
			dProton_FDCdEdx_P->Fill(locProtonP4.Vect().Mag(), FDC_dEdx);
			dProton_P_Theta->Fill(locProtonP4.Vect().Theta()*180/TMath::Pi(), locProtonP4.Vect().Mag());

			dPion_CDCdEdx_P->Fill(locPion1P4.Vect().Mag(), CDC_dEdx_pi1);
			dPion_FDCdEdx_P->Fill(locPion1P4.Vect().Mag(), FDC_dEdx_pi1);

            double locDeltaPhi = (locProtonP4.Phi() - locPhiP4.Phi())*180./TMath::Pi();
            if(locDeltaPhi > 360.) locDeltaPhi -= 360.;
            if(locDeltaPhi < 0.) locDeltaPhi += 360.;
            dDeltaPhi->Fill(locDeltaPhi);

			if(fabs(locP4_2pi.M() - kK0Mass) < 0.1) {
				dMKSKL->Fill(locPhiP4.M());

                vector<const DMCThrown *> mcthrowns;
                locEventLoop->Get(mcthrowns);
                //locEventLoop->Get(mcthrowns, "Primary");
                if(mcthrowns.size() > 0) {
                    cout << "Event!" << endl;
                    for(int loc_i=0; loc_i<mcthrowns.size(); loc_i++) {
                        vector<pair<string,string> >  items;                
                        mcthrowns[loc_i]->toStrings(items);
                        cout << "#" << (loc_i+1) << " ->   ";
                        for(int j=0; j<items.size(); j++) 
                            cout << items[j].first << " : " << items[j].second << ", ";
                        cout << endl;
                    }
                        
                }

                vector<const DNeutralParticle *> photons;
                locEventLoop->Get(photons);
                if(photons.size()>=2) {
                    for(int loc_gam_i=0; loc_gam_i<photons.size()-1; loc_gam_i++) {
                        DLorentzVector gam1 = photons[loc_gam_i]->Get_Hypothesis(Gamma)->lorentzMomentum();
                        for(int loc_gam_j=loc_gam_i+1; loc_gam_j<photons.size(); loc_gam_j++) {
                            DLorentzVector gam2 = photons[loc_gam_j]->Get_Hypothesis(Gamma)->lorentzMomentum();
                            dMgg->Fill( (gam1+gam2).M() );
                        }
                    }
                }
			}
			if(fabs(locP4_2pi.M() - kK0Mass) < 0.05) {
				dMKSKL2->Fill(locPhiP4.M());
			}
			if(fabs(locP4_2pi.M() - kK0Mass) < 0.02) {
				dMKSKL3->Fill(locPhiP4.M());
			}

		}

        // ---
		if( (fabs(locPhiP4.M() - kPhiMass) < 0.07) 
		    && (fabs(locP4_2pi.M() - kK0Mass)  < 0.1) ) {

			if(fabs(locP4_2pi.M() - kK0Mass) < 0.1) {
				dMKL->Fill(locMissingP4.M());
			}
			if(fabs(locP4_2pi.M() - kK0Mass) < 0.05) {
				dMKL2->Fill(locMissingP4.M());
			}
			if(fabs(locP4_2pi.M() - kK0Mass) < 0.02) {
				dMKL3->Fill(locMissingP4.M());
			}
            
        }
	}
	Unlock_Action(); //RELEASE ROOT LOCK!!
    */


	//Minor change to push to svn

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