#include "DAnalysisActionEtaPrime.h"
#include <vector>
#include <GlueX.h>
#include <deque>
#include <PID/DKinematicData.h>

void DAnalysisActionEtaPrime::Initialize(JEventLoop* locEventLoop)
{
  // JApplication *japp=locEventLoop->GetJApplication();

  myevt=0;

  //CREATE THE HISTOGRAMS
  // Get_Application()->RootWriteLock(); //ACQUIRE ROOT LOCK!!
  japp->RootWriteLock();
  CreateAndChangeTo_ActionDirectory();
  
  h_p_missing_mass=static_cast<TH1D*>(gDirectory->Get("h_p_missing_mass"));
  if(!h_p_missing_mass){
    h_p_missing_mass=new TH1D("h_p_missing_mass","Missing mass off proton",
			      250, 0., 2.5);
    h_p_missing_mass->SetXTitle("Missing mass [GeV]");
  } 

  h_theta_vs_p=static_cast<TH2D*>(gDirectory->Get("h_theta_vs_p"));
  if(!h_theta_vs_p){
    h_theta_vs_p=new TH2D("h_theta_vs_p","Proton #theta_{LAB} vs. p",
			  200,0,2.,180,0.,90.);
    h_theta_vs_p->SetXTitle("p [GeV/c]");
    h_theta_vs_p->SetYTitle("#theta [degrees]");
  }

  h_2gamma_2d=static_cast<TH2D*>(gDirectory->Get("h_2gamma_2d"));
  if (!h_2gamma_2d){
    h_2gamma_2d=new TH2D("h_2gamma_2d","2#gamma pair 2 mass vs pair 1 mass",100,0,1,100,0,1);
    h_2gamma_2d->SetXTitle("m(2#gamma) pair 1 [Gev]");
    h_2gamma_2d->SetYTitle("m(2#gamma) pair 2 [Gev]");
  }

  gDirectory->cd("..");

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


bool DAnalysisActionEtaPrime::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
{
      
  // Results from kinematic fit
  const DParticleComboStep *comboStep=locParticleCombo->Get_ParticleComboStep(0);
  DLorentzVector proton=comboStep->Get_FinalParticle(0)->lorentzMomentum();
  DLorentzVector beam = comboStep->Get_InitialParticle()->lorentzMomentum();
  DLorentzVector target =comboStep->Get_TargetParticle()->lorentzMomentum();
  
  DLorentzVector eta=beam+target-proton;
  double eta_mass=eta.M(); 

  deque<const DKinematicData*>pims;
  comboStep->Get_FinalParticles(PiMinus,pims);
  deque<const DKinematicData*>pips;
  comboStep->Get_FinalParticles(PiPlus,pips);
     
  deque<const DKinematicData*>neutrals;
  comboStep->Get_FinalParticles(Gamma,neutrals);
  
  deque<const DKinematicData*>measured_neutrals;
  comboStep->Get_FinalParticles_Measured(Gamma,measured_neutrals);
  
  vector<DLorentzVector>measured_gams(neutrals.size());
  for (unsigned int i=0;i<measured_neutrals.size();i++){
    measured_gams[i]=measured_neutrals[i]->lorentzMomentum();
  }
  
  
  vector<DLorentzVector>gams(neutrals.size());
  for (unsigned int i=0;i<neutrals.size();i++){
    gams[i]=neutrals[i]->lorentzMomentum();
  }
  
  Lock_Action(); //ACQUIRE ROOT LOCK!!
  {
    //      h_p_missing_mass->Fill(rareEtaCombos[0].eta_mass);
    h_theta_vs_p->Fill(proton.P(),(180./M_PI)*proton.Theta());
    h_p_missing_mass->Fill(eta_mass);
    
    if (gams.size()==4){
      DLorentzVector pair1=gams[0]+gams[1];
      DLorentzVector pair2=gams[2]+gams[3];
      h_2gamma_2d->Fill(pair1.M(),pair2.M());
      pair1=gams[0]+gams[2];
      pair2=gams[1]+gams[3];
      h_2gamma_2d->Fill(pair1.M(),pair2.M());
      pair1=gams[0]+gams[3];
      pair2=gams[1]+gams[2];
      h_2gamma_2d->Fill(pair1.M(),pair2.M());
    }
  }
  Unlock_Action(); //RELEASE ROOT LOCK!!


  return true;
}