#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_t_pi0=static_cast<TH1D*>(gDirectory->Get("h_t_pi0"));
  if(!h_t_pi0){
    h_t_pi0=new TH1D("h_t_pi0","#pi0 t distribution",
			      1000, 0., 1.0);
    h_t_pi0->SetXTitle("-t [GeV^2]");
  } 
  h_t_eta=static_cast<TH1D*>(gDirectory->Get("h_t_eta"));
  if(!h_t_eta){
    h_t_eta=new TH1D("h_t_eta","#eta t distribution",
			      1000, 0., 1.0);
    h_t_eta->SetXTitle("-t [GeV^2]");
  } 
  h_p_eta_mass=static_cast<TH1D*>(gDirectory->Get("h_p_eta_mass"));
  if(!h_p_eta_mass){
    h_p_eta_mass=new TH1D("h_p_eta_mass","proton#eta mass",
			      200, 0., 4.0);
    h_p_eta_mass->SetXTitle("Mass [GeV]");
  } 
  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",
			      200, 0., 2.0);
    h_p_missing_mass->SetXTitle("Missing mass [GeV]");
  } 
  h_p_missing_mass_cut=static_cast<TH1D*>(gDirectory->Get("h_p_missing_mass_cut"));
  if(!h_p_missing_mass_cut){
    h_p_missing_mass_cut=new TH1D("h_p_missing_mass_cut","Missing mass off proton, #pi0 required",
			      200, 0., 2.0);
    h_p_missing_mass_cut->SetXTitle("Missing mass [GeV]");
  }
   h_p_missing_mass_eta_cut=static_cast<TH1D*>(gDirectory->Get("h_p_missing_mass_eta_cut"));
   if(!h_p_missing_mass_eta_cut){
    h_p_missing_mass_eta_cut=new TH1D("h_p_missing_mass_eta_cut","Missing mass off proton, #eta required",
			      200, 0., 2.0);
    h_p_missing_mass_eta_cut->SetXTitle("Missing mass [GeV]");
  }
  h_2gamma_invariant_mass=static_cast<TH1D*>(gDirectory->Get("h_2gamma_invariant_mass"));
  if(!h_2gamma_invariant_mass){
    h_2gamma_invariant_mass=new TH1D("h_2gamma_invariant_mass","2#gamma invariant mass",
				     200, 0.0, 2.0);
    h_2gamma_invariant_mass->SetXTitle("M(2#gamma) [GeV]");
  }
  h_2gamma_measured_invariant_mass=static_cast<TH1D*>(gDirectory->Get("h_2gamma_measured_invariant_mass"));
  if(!h_2gamma_measured_invariant_mass){
    h_2gamma_measured_invariant_mass=new TH1D("h_2gamma_measured_invariant_mass","2#gamma invariant mass",
				     200, 0.0, 2.0);
    h_2gamma_measured_invariant_mass->SetXTitle("M(2#gamma) [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_p_mm_vs_E=static_cast<TH2D*>(gDirectory->Get("h_p_mm_vs_E"));
  if(!h_p_mm_vs_E){
    h_p_mm_vs_E=new TH2D("h_p_mm_vs_E","Missing mass off proton  vs. E#gamma",
			  100,0,6.,100,0.,1.5);
    h_p_mm_vs_E->SetXTitle("E#gamma [GeV]");
    h_p_mm_vs_E->SetYTitle("Missing mass [GeV]");
  }
  h_2gamma_dphi=static_cast<TH1D*>(gDirectory->Get("h_2gamma_dphi"));
  if (!h_2gamma_dphi){
    h_2gamma_dphi=new TH1D("h_2gamma_dphi","#delta#phi (2#gamma)",200,-6.5,6.5);
    h_2gamma_dphi->SetTitle("#delta#phi [rad]");
  }
  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]");
  }
  h_pip_pim_mass=static_cast<TH1D*>(gDirectory->Get("h_pip_pim_mass"));
  if (!h_pip_pim_mass){
    h_pip_pim_mass=new TH1D("h_pip_pim_mass","#pi+#pi- mass",200,0,2);
    h_pip_pim_mass->SetTitle("M(#pi+#pi-) [GeV]");
  }
  h_ppippim_missing_mass_sq=static_cast<TH1D*>(gDirectory->Get("h_ppippim_missing_mass_sq"));
  if (!h_ppippim_missing_mass_sq){
    h_ppippim_missing_mass_sq=new TH1D("h_ppippim_missing_mass_sq","p#pi+#pi- missing mass squared",200,-0.5,1.5);
    h_ppippim_missing_mass_sq->SetTitle("MM^{2}(p#pi+#pi-) [GeV^{2}]");
  }
  gDirectory->cd("..");

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


bool DAnalysisActionEtaPrime::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
{
  
  unsigned int evtno=locEventLoop->GetJEvent().GetEventNumber();
      
  // 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(); 
  DVector3 boost=-eta.BoostVector();
  double t=(beam-eta).M2();

  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();
  }
  
  japp->RootWriteLock();
  
  if (pips.size()==1 && pims.size()==1){
    DLorentzVector pip_pim=pips[0]->lorentzMomentum()+pims[0]->lorentzMomentum();
    h_pip_pim_mass->Fill(pip_pim.M());
    DLorentzVector ppippim_missing=eta-pip_pim;
    h_ppippim_missing_mass_sq->Fill(ppippim_missing.M2());
  }
  
  for (unsigned int i=0;i<measured_gams.size();i++){
    for (unsigned int j=i+1;j<measured_gams.size();j++){
      DLorentzVector twogam=measured_gams[i]+measured_gams[j];
      
      double twogamma_mass=twogam.M();
      
      h_2gamma_measured_invariant_mass->Fill(twogamma_mass);
    }
  }
  
  //      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);
  bool got_pi0=false;
  bool got_eta=false;
  for (unsigned int i=0;i<gams.size();i++){
    for (unsigned int j=i+1;j<gams.size();j++){
      DLorentzVector twogam=gams[i]+gams[j];
      
      double twogamma_mass=twogam.M();
      gams[i].Boost(boost);
      gams[j].Boost(boost);
      double dphi_rest=gams[i].Phi()-gams[j].Phi();
      if (dphi_rest<-M_PI)dphi_rest+=2*M_PI;
      if (dphi_rest>M_PI) dphi_rest-=2*M_PI;
      h_2gamma_dphi->Fill(dphi_rest);
      
      
      h_2gamma_invariant_mass->Fill(twogamma_mass);
      if (twogamma_mass>0.1&&twogamma_mass<0.18) got_pi0=true;
      if (twogamma_mass>0.5&&twogamma_mass<0.6) got_eta=true;
    }
  }
  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());
  }
  
  if (got_pi0){
    h_p_missing_mass_cut->Fill(eta_mass);
    h_t_pi0->Fill(-t);
  }
  if (got_eta){
    h_p_missing_mass_eta_cut->Fill(eta_mass);
    h_p_eta_mass->Fill((beam+target).M());
    h_t_eta->Fill(-t);
  }
  h_p_mm_vs_E->Fill(beam.E(),eta_mass);
  

  japp->RootUnLock();

  return true;
}