#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);
  }
  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)
{
      // 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;
      comboStep=locParticleCombo->Get_ParticleComboStep(1);
      for (unsigned int i=0;i<3;i++){
	pi[i]= comboStep->Get_FinalParticle(i)->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;
  
      Lock_Action(); //ACQUIRE ROOT LOCK!!
      {
	if (eta_mass>0.53 && eta_mass<0.58){
	  dalitzXY->Fill(x_dalitz,y_dalitz);
	  dalitzZ->Fill(z_dalitz);
	}
      }
      Unlock_Action(); //RELEASE ROOT LOCK!!

      

  return true;
}