#include "DAnalysisActionEta.h"
#include <FCAL/DFCALShower.h>
#include <vector>

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

  myevt=0;

  //CREATE THE HISTOGRAMS
  // Get_Application()->RootWriteLock(); //ACQUIRE ROOT LOCK!!
  japp->RootWriteLock();
  CreateAndChangeTo_ActionDirectory();

  Htwo_gamma_mass = static_cast<TH1D*>(gDirectory->Get("Htwo_gamma_mass"));
  if (!Htwo_gamma_mass){
    Htwo_gamma_mass=new TH1D("Htwo_gamma_mass","Invariant Mass",100,0.,1.0);
    Htwo_gamma_mass->SetXTitle("M(#gamma#gamma) [GeV]");
    Htwo_gamma_mass->SetYTitle("counts/0.005 GeV");
  } 
  Hpi0_two_gamma_mass = static_cast<TH1D*>(gDirectory->Get("Hpi0_two_gamma_mass"));
  if (!Hpi0_two_gamma_mass){
    Hpi0_two_gamma_mass=new TH1D("Hpi0_two_gamma_mass","Invariant Mass",100,0.,2.0);
    Hpi0_two_gamma_mass->SetXTitle("M(#pi0#gamma#gamma) [GeV]");
    Hpi0_two_gamma_mass->SetYTitle("counts/0.02 GeV");
  }
  Hpi0_pi0_mass = static_cast<TH1D*>(gDirectory->Get("Hpi0_pi0_mass"));
  if (!Hpi0_pi0_mass){
    Hpi0_pi0_mass=new TH1D("Hpi0_pi0_mass","Invariant Mass",100,0.,2.0);
    Hpi0_pi0_mass->SetXTitle("M(#pi0#pi0) [GeV]");
    Hpi0_pi0_mass->SetYTitle("counts/0.02 GeV");
  }
  Heta_pi0_mass = static_cast<TH1D*>(gDirectory->Get("Heta_pi0_mass"));
  if (!Heta_pi0_mass){
    Heta_pi0_mass=new TH1D("Heta_pi0_mass","Invariant Mass",100,0.,2.0);
    Heta_pi0_mass->SetXTitle("M(#eta#pi0) [GeV]");
    Heta_pi0_mass->SetYTitle("counts/0.02 GeV");
  }
  Hpi0g_dalitz = static_cast<TH2D*>(gDirectory->Get("Hpi0g_dalitz"));
  if (!Hpi0g_dalitz){
    Hpi0g_dalitz=new TH2D("Hpi0g_dalitz","#pi^{0}#gamma Dalitz plot",100,0.,0.35,100,0,0.35);
    Hpi0g_dalitz->SetXTitle("M^{2}(#pi^{0}#gamma_{1}) [GeV^{2}]");
    Hpi0g_dalitz->SetYTitle("M^{2}(#pi^{0}#gamma_{2}) [GeV^{2}]");
  }

  gDirectory->cd("..");

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

bool DAnalysisActionEta::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
{
      // Results from kinematic fit
      DLorentzVector proton=locParticleCombo->Get_ParticleComboStep(0)->Get_FinalParticle(0)->lorentzMomentum();
      DLorentzVector beam = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle()->lorentzMomentum();
      DLorentzVector target =locParticleCombo->Get_ParticleComboStep(0)->Get_TargetParticle()->lorentzMomentum();
      
      DLorentzVector pi0 = locParticleCombo->Get_ParticleComboStep(0)->Get_FinalParticle(1)->lorentzMomentum();
      DLorentzVector gamma1 = locParticleCombo->Get_ParticleComboStep(0)->Get_FinalParticle(2)->lorentzMomentum();
      DLorentzVector gamma2 = locParticleCombo->Get_ParticleComboStep(0)->Get_FinalParticle(3)->lorentzMomentum();
      DLorentzVector two_gamma=gamma1+gamma2;
      DLorentzVector pi0_g1=pi0+gamma1;
      DLorentzVector pi0_g2=pi0+gamma2;
      DLorentzVector pi0_2gamma=pi0_g1+gamma2;
      double two_gamma_mass=two_gamma.M();
      double four_gamma_mass=pi0_2gamma.M();
      
      japp->RootWriteLock();
       
      Htwo_gamma_mass->Fill(two_gamma_mass);
      
      if (two_gamma_mass>0.1&&two_gamma_mass<0.18){
	Hpi0_pi0_mass->Fill(four_gamma_mass);
      }
      if (two_gamma_mass>0.5&&two_gamma_mass<0.6){
	Heta_pi0_mass->Fill(four_gamma_mass);
      }
      Hpi0g_dalitz->Fill(pi0_g1.M2(),pi0_g2.M2());
      Hpi0_two_gamma_mass->Fill(four_gamma_mass);
      
      japp->RootUnLock();
    
  
    
  //myevt=evtno;

  return true;
}