#include "DAnalysisActionEtaDalitz.h"
#include <FCAL/DFCALShower.h>
#include <ANALYSIS/DParticleComboStep.h>
#include <vector>

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

  myevt=0;

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

  dalitzZ = static_cast<TH1D*>(gDirectory->Get("dalitzZ"));
  if (!dalitzZ){
    dalitzZ=new TH1D("dalitzZ","dalitz Z",110,-0.05,1.05);
  }
  dalitzXY=static_cast<TH2D*>(gDirectory->Get("dalitzXY"));
  if (!dalitzXY){
    dalitzXY=new TH2D("dalitzXY","Dalitz distribution Y vs X",100,-1.,1.,100,-1.,1);
  }
  hEgamma=static_cast<TH1D*>(gDirectory->Get("hEgamma"));
  if (!hEgamma){
    hEgamma=new TH1D("hEgamma","E_{#gamma} distribution",1000,0,12.);
    hEgamma->SetTitle("E_{#gamma} [GeV]");
  } 
  hInvMass=static_cast<TH1D*>(gDirectory->Get("hInvMass"));
  if (!hInvMass){
    hInvMass=new TH1D("hInvMass","3 pion mass",200,0,2.);
    hInvMass->SetTitle("m{3#pi} [GeV]");
  }
  htheta_vs_p=static_cast<TH2D*>(gDirectory->Get("htheta_vs_p"));
  if(!htheta_vs_p){
    htheta_vs_p=new TH2D("htheta_vs_p","Proton #theta_{LAB} vs. p",
			       200,0,2.,180,0.,90.);
    htheta_vs_p->SetXTitle("p [GeV/c]");
    htheta_vs_p->SetYTitle("#theta [degrees]");
  }
			     
  gDirectory->cd("..");

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


bool dalitzCombo_cmp(DAnalysisActionEtaDalitz::DalitzInfo a,
		     DAnalysisActionEtaDalitz::DalitzInfo b){
  return a.cl>b.cl;
}

bool DAnalysisActionEtaDalitz::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 pi[3];
      DLorentzVector eta;
      for (unsigned int i=0;i<3;i++){
	pi[i]= comboStep->Get_FinalParticle(i+1)->lorentzMomentum();
	eta+=pi[i];
      }
      
      DVector3 boost=-eta.BoostVector();
      double eta_mass=eta.M();
      // Boost to the eta rest frame
      for (unsigned int i=0;i<3;i++){
	pi[i].Boost(boost);
      }
    
      double m1=pi[0].M(),m2=pi[1].M(),m3=pi[2].M();
      double T1=pi[0].E()-m1; // pi0 for charged channel
      double T2=pi[1].E()-m2; // pi+ for charged channel
      double T3=pi[2].E()-m3; // pi- for charged channel
      double Q_eta=eta_mass-m1-m2-m3;
      double x_dalitz=sqrt(3.)*(T2-T3)/Q_eta;
      double y_dalitz=3.*T1/Q_eta-1.;
      double z_dalitz=x_dalitz*x_dalitz+y_dalitz*y_dalitz;
  
      japp->RootWriteLock();
      
      if (eta_mass>0.53 && eta_mass<0.58){
	dalitzXY->Fill(x_dalitz,y_dalitz);
	dalitzZ->Fill(z_dalitz);
      }
      hEgamma->Fill(beam.E());
      hInvMass->Fill(eta_mass);
      htheta_vs_p->Fill(proton.P(),(180./M_PI)*proton.Theta());
      
      japp->RootUnLock();

  myevt=evtno;

  return true;
}