#include "DAnalysisActionEtaU.h"
#include <FCAL/DFCALShower.h>

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

  DApplication* locApplication = dynamic_cast<DApplication*>(japp);
  //  dKinFitter.Set_BField(locApplication->GetBfieldPointer());


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

  uboson_mass = static_cast<TH1D*>(gDirectory->Get("uboson_mass"));
  if (!uboson_mass){
    uboson_mass=new TH1D("uboson_mass","e+e- mass",1000,0.,1.0);
    uboson_mass->SetXTitle("M(e^{+}e^{-}) [GeV]");
    uboson_mass->SetYTitle("counts/1 MeV");
  }
  uboson_mass_kin = static_cast<TH1D*>(gDirectory->Get("uboson_mass_kin"));
  if (!uboson_mass_kin){
    uboson_mass_kin=new TH1D("uboson_mass_kin","e+e- mass, momenta from kinematic fit",1000,0.,1.0);
    uboson_mass_kin->SetXTitle("M(e^{+}e^{-}) [GeV]");
    uboson_mass_kin->SetYTitle("counts/1 MeV");
  }


  p_missing_mass = static_cast<TH1D*>(gDirectory->Get("p_missing_mass"));
  if (!p_missing_mass){
    p_missing_mass=new TH1D("p_missing_mass","missing mass off proton",1000,0.,1.0);
    p_missing_mass->SetXTitle("MM(#gammap->p(X)) [GeV]");
    p_missing_mass->SetYTitle("counts/1 MeV");
  }

  ep_em_gamma_mass = static_cast<TH1D*>(gDirectory->Get("ep_em_gamma_mass"));
  if (!ep_em_gamma_mass){
    ep_em_gamma_mass=new TH1D("ep_em_gamma_mass","e^{+}e^{-}#gamma mass",1000,0.,1.0);
    ep_em_gamma_mass->SetXTitle("M(e^{+}e^{-}#gamma) [GeV]");
    ep_em_gamma_mass->SetYTitle("counts/1 MeV");
  }

  ep_em_gamma_p_MM2 = static_cast<TH1D*>(gDirectory->Get("ep_em_gamma_p_MM2"));
  if (!ep_em_gamma_p_MM2){
    ep_em_gamma_p_MM2=new TH1D("ep_em_gamma_p_mm2","e^{+}e^{-}#gammap missing mass squared.",1000,-0.5,0.5); 
    ep_em_gamma_p_MM2->SetXTitle("MM^{2}(e^{+}e^{-}#gammap) [GeV^{2}]");
    ep_em_gamma_p_MM2->SetYTitle("counts/0.001 GeV^{2}");
  }

   ep_em_gamma_p_missing_E = static_cast<TH1D*>(gDirectory->Get("ep_em_gamma_p_missing_E"));
  if (!ep_em_gamma_p_missing_E){
    ep_em_gamma_p_missing_E=new TH1D("ep_em_gamma_p_missing_E","e^{+}e^{-}#gammap missing energy",200,-0.5,9.5); 
    ep_em_gamma_p_missing_E->SetXTitle("Missing Energy [GeV]");
    ep_em_gamma_p_missing_E->SetYTitle("counts");
  }

  uboson_mass_cut = static_cast<TH1D*>(gDirectory->Get("uboson_mass_cut"));
  if(!uboson_mass_cut){
    uboson_mass_cut=new TH1D("uboson_mass_cut","e+e- mass, e^{+}e^{-} E/p  cut",1000,0.,1.0);  
    uboson_mass_cut->SetXTitle("M(e^{+}e^{-}) [GeV]");
    uboson_mass_cut->SetYTitle("counts/1 MeV");
  }

  ep_em_gamma_mass_cut = static_cast<TH1D*>(gDirectory->Get("ep_em_gamma_mass_cut"));
  if(!ep_em_gamma_mass_cut){
    ep_em_gamma_mass_cut=new TH1D("ep_em_gamma_mass_cut","e^{+}e^{-}#gamma mass, e^{+}e^{-} E/p cut",1000,0.,1.0);
    ep_em_gamma_mass_cut->SetXTitle("M(e^{+}e^{-}#gamma) [GeV]");
    ep_em_gamma_mass_cut->SetYTitle("counts/1 MeV");
  }
  ep_em_gamma_mass_kin = static_cast<TH1D*>(gDirectory->Get("ep_em_gamma_mass_kin"));
  if(!ep_em_gamma_mass_kin){
    ep_em_gamma_mass_kin=new TH1D("ep_em_gamma_mass_kin","e^{+}e^{-}#gamma mass, e^{+}e^{-} E/p cut, kinematic fit",1000,0.,1.0);
    ep_em_gamma_mass_kin->SetXTitle("M(e^{+}e^{-}#gamma) [GeV]");
    ep_em_gamma_mass_kin->SetYTitle("counts/1 MeV");
  }

  ep_em_gamma_p_MM2_cut = static_cast<TH1D*>(gDirectory->Get("ep_em_gamma_p_MM2_cut"));
  if (!ep_em_gamma_p_MM2_cut){
    ep_em_gamma_p_MM2_cut=new TH1D("ep_em_gamma_p_mm2_cut","e^{+}e^{-}#gammap missing mass squared, e^{+}e^{-} E/p cut",1000,-0.5,0.5); 
    ep_em_gamma_p_MM2_cut->SetXTitle("MM^{2}(e^{+}e^{-}#gammap) [GeV^{2}]");
    ep_em_gamma_p_MM2_cut->SetYTitle("counts/0.001 GeV^{2}");
  }

  mandelstam_t = static_cast<TH1D*>(gDirectory->Get("mandelstam_t"));
  if (!mandelstam_t){
    mandelstam_t=new TH1D("mandelstam_t","t",100,-2.,0.);
    mandelstam_t->SetXTitle("t [GeV^{2}]");
    mandelstam_t->SetYTitle("counts/0.02 GeV^{2}");
  }
  mandelstam_t_kin = static_cast<TH1D*>(gDirectory->Get("mandelstam_t_kin"));
  if (!mandelstam_t_kin){
    mandelstam_t_kin=new TH1D("mandelstam_t_kin","t, after kinematic fit",100,-2.,0.);
    mandelstam_t_kin->SetXTitle("t [GeV^{2}]");
    mandelstam_t_kin->SetYTitle("counts/0.02 GeV^{2}");
  }

  electron_E_over_p=static_cast<TH2D*>(gDirectory->Get("electron_E_over_p"));
  if (!electron_E_over_p){
    electron_E_over_p=new TH2D("electron_E_over_p","E/p vs p for e^{-}",100,0,10,120,0.0,1.2);
    electron_E_over_p->SetXTitle("p [GeV/c]");
    electron_E_over_p->SetYTitle("E/p");
  }
	 
  positron_E_over_p=static_cast<TH2D*>(gDirectory->Get("positron_E_over_p"));
  if (!positron_E_over_p){
    positron_E_over_p=new TH2D("positron_E_over_p","E/p vs p for e^{+}",100,0,10,120,0.0,1.2);			       
    positron_E_over_p->SetXTitle("p [GeV/c]");
    positron_E_over_p->SetYTitle("E/p");
  }

  ep_em_vertex_hist=static_cast<TH2D*>(gDirectory->Get("ep_em_vertex_hist"));
  if (!ep_em_vertex_hist){
    ep_em_vertex_hist=new TH2D("ep_em_vertex_hist","e^{+}e^{-} vertex r vs z",500,0,100,200,0,10);
    ep_em_vertex_hist->SetXTitle("z [cm]");
    ep_em_vertex_hist->SetYTitle("r [cm]");
  }
  vertex_hist=static_cast<TH2D*>(gDirectory->Get("vertex_hist"));
  if (!vertex_hist){
    vertex_hist=new TH2D("vertex_hist","e^{+}e^{-}p vertex r vs z",500,0,100,200,0,10);
    vertex_hist->SetXTitle("z [cm]");
    vertex_hist->SetYTitle("r [cm]");
  }		
  
  ep_em_vertex_kin=static_cast<TH2D*>(gDirectory->Get("ep_em_vertex_kin"));
  if (!ep_em_vertex_kin){
    ep_em_vertex_kin=new TH2D("ep_em_vertex_kin","e^{+}e^{-} kinematic fit vertex r vs z",500,0,100,200,0,10);
    ep_em_vertex_kin->SetXTitle("z [cm]");
    ep_em_vertex_kin->SetYTitle("r [cm]");
  }
  vertex_kin=static_cast<TH2D*>(gDirectory->Get("vertex_kin"));
  if (!vertex_kin){
    vertex_kin=new TH2D("vertex_kin","e^{+}e^{-}p kinematic fit vertex r vs z",500,0,100,200,0,10);
    vertex_kin->SetXTitle("z [cm]");
    vertex_kin->SetYTitle("r [cm]");
  }		

  gDirectory->cd("..");

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

}

bool DAnalysisActionEtaU::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo)
{
     const DChargedTrackHypothesis* electron = static_cast<const DChargedTrackHypothesis*>(locParticleCombo->Get_ParticleComboStep(1)->Get_FinalParticle(0));
    const DChargedTrackHypothesis* positron = static_cast<const DChargedTrackHypothesis*>(locParticleCombo->Get_ParticleComboStep(1)->Get_FinalParticle(1));
    const DChargedTrackHypothesis* proton = static_cast<const DChargedTrackHypothesis*>(locParticleCombo->Get_ParticleComboStep(0)->Get_FinalParticle(1));
    const DNeutralParticleHypothesis* photon = static_cast<const DNeutralParticleHypothesis*>(locParticleCombo->Get_ParticleComboStep(1)->Get_FinalParticle(2));
    const DBeamPhoton* beam = static_cast<const DBeamPhoton*>(locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle());
    
    DLorentzVector beam_p4 = beam->lorentzMomentum();
    DLorentzVector photon_p4=photon->lorentzMomentum();
    DLorentzVector target_p4(0.,0.,0.,0.93827);
    DLorentzVector electron_p4= electron->lorentzMomentum();
    DLorentzVector positron_p4= positron->lorentzMomentum(); 
    DLorentzVector proton_p4= proton->lorentzMomentum();
    DLorentzVector uboson_p4=electron_p4+positron_p4;
    
    DLorentzVector init_reaction_p4=beam_p4+target_p4;
    DLorentzVector proton_missing_p4=init_reaction_p4-proton_p4;
    DLorentzVector eta_p4=electron_p4+positron_p4+photon_p4;
    DLorentzVector eta_p_missing_p4=proton_missing_p4-eta_p4;
    double t=(beam_p4-eta_p4).M2();
    
    vector<const DFCALShower*>matched_electron_showers;
    electron->GetT(matched_electron_showers); 
    vector<const DFCALShower*>matched_positron_showers;
    positron->GetT(matched_positron_showers);
    
    DVector3 boost=-eta_p4.BoostVector();
    DLorentzVector temp=eta_p4;
    temp.Boost(boost);
    
    //  JApplication *japp=locEventLoop->GetJApplication();
    japp->RootWriteLock();
    
    bool good_electron=false;
    bool good_positron=false;
    if (matched_electron_showers.size()){
      double p=electron_p4.P();
      double E_over_p=matched_electron_showers[0]->getEnergy()/p;
      double one_minus_E_over_p=1.-E_over_p;
      double sig_electron=0.0646/sqrt(p)+0.0109*p-0.0172;
      good_electron=(fabs(one_minus_E_over_p)<3.*sig_electron)?true:false;
      electron_E_over_p->Fill(p,E_over_p);
    }
    
    if (matched_positron_showers.size()){
      double p=positron_p4.P();
      double E_over_p=matched_positron_showers[0]->getEnergy()/p;  
      double one_minus_E_over_p=1.-E_over_p;
      double sig_positron=0.0675/sqrt(p)+0.0116*p-0.0204;
      good_positron=(fabs(one_minus_E_over_p)<3.*sig_positron)?true:false;
      positron_E_over_p->Fill(p,E_over_p);
    }
    
    if (good_electron && good_positron /* && eta_p_missing_p4.E()<1.0*/){
      // Histograms using reconstructed 4-momenta, pre-kinematic fit
      // Missing energy
      ep_em_gamma_p_missing_E->Fill(eta_p_missing_p4.E());
      
      if (eta_p_missing_p4.E()<1.0){
	mandelstam_t->Fill(t);
	uboson_mass->Fill(uboson_p4.M());
	ep_em_gamma_mass->Fill(eta_p4.M());
	ep_em_gamma_p_MM2->Fill(eta_p_missing_p4.M2());
	/*
	DKinematicData electron_kd,positron_kd;
	DVector3 ep_em_vertex;
	double doca,var_doca;
	electron->dRT->IntersectTracks(positron->dRT,&electron_kd,&positron_kd,ep_em_vertex,doca,var_doca);
	ep_em_vertex_hist->Fill(ep_em_vertex.z(),ep_em_vertex.Perp());
	
	DKinematicData proton_kd;
	proton->dRT->FindPOCAtoPoint(ep_em_vertex,NULL,&proton_kd,doca,var_doca);
	DVector3 vertex=0.5*(ep_em_vertex+proton_kd.position());
	vertex_hist->Fill(vertex.z(),vertex.Perp());
	*/
	p_missing_mass->Fill(proton_missing_p4.M());
	
      }
    }

  
  japp->RootUnLock();
 
  return true;
}