#include <iostream>
using namespace std;

#include "MyProcessor.h"
#include "HDDM/hddm_s.h"
#include "HDDM/DEventSourceHDDM.h"

#include "TRACKING/DMCThrown.h"
#include "PID/DPhoton.h"
//#include "FCAL/DFCALPhoton.h"
#include "PID/DBeamPhoton.h"
//#include "PID/DTwoGammaFit_factory_PI0.h"
//#include "PID/DParticle.h"

#include <TDirectory.h>
#include <TDatabasePDG.h>


void MyProcessor::center_axis(TH1* hist)
{
  hist->GetXaxis()->CenterTitle(kTRUE);
  hist->GetYaxis()->CenterTitle(kTRUE);
}

//=========================================
void MyProcessor::center_axis(TH1* hist, TDirectory* dir)
{
  hist->SetDirectory(dir);
  hist->GetXaxis()->CenterTitle(kTRUE);
  hist->GetYaxis()->CenterTitle(kTRUE);
}

//============================================

void MyProcessor::center_axis(TH2* hist)
{
  hist->GetXaxis()->CenterTitle(kTRUE);
  hist->GetYaxis()->CenterTitle(kTRUE);
}

//===========================================
void MyProcessor::center_axis(TH2* hist, TDirectory* dir)
{
  hist->SetDirectory(dir);
  hist->GetXaxis()->CenterTitle(kTRUE);
  hist->GetYaxis()->CenterTitle(kTRUE);
}

//------------------------------------------------------------------
// init   -Open output file here (e.g. a ROOT file)
//------------------------------------------------------------------
jerror_t MyProcessor::init(void)
{
  // open ROOT file
  ROOTfile = new TFile("hd_ana.root","RECREATE","Produced by hd_ana");
  cout<<"Opened ROOT file \"hd_ana.root\""<<endl;
  
  //->>>General Histograms 
  //-------------------------
  E_beam = new TH1F("E_beam", "; E_{#gamma} (GeV); counts/1.0MeV", 12000, 0., 12.);
  center_axis(E_beam);

  fcal_time = new TH1F("fcal_time", "; FCAL time (ns); counts/0.5ns", 480, -120., 120.);
  center_axis(fcal_time);

  xy_position = new TH2F("xy_position", "; X_{FCAL} (cm); Y_{FCAL} (cm)", 2000, -100., 100., 2000, -100., 100.);
  center_axis(xy_position);

  num_photons = new TH1F("num_photons", "; number of #gamma's in FCAL per event; counts/event", 14, -0.5, 13.5);
  center_axis(num_photons);

  num_charged = new TH1F("num_charged", "; number of  charged tracks per event; counts/event", 14, -0.5, 13.5);
  center_axis(num_charged);

  total_energy_vs_elast = new TH2F("total_energy_vs_elast", "; E_{total, FCAL}/E_{#gamma}; E_{total, FCAL} (GeV)", 300, 0., 1.5, 12000, 0., 12.);
  center_axis(total_energy_vs_elast);

  //  //->>>>>>>General Histograms N_FCAL_PHOTONS>2; E_FCAL_Total>9.0
  //
  fcal_time_cut = new TH1F("fcal_time_cut", "; FCAL time (ns); counts/0.5ns", 480, -120., 120.);
  center_axis(fcal_time_cut);

  xy_position_cut = new TH2F("xy_position_cut", "; X_{FCAL} (cm); Y_{FCAL} (cm)", 2000, -100., 100., 2000, -100., 100.);
  center_axis(xy_position);

  num_photons_cut = new TH1F("num_photons_cut", "; number of #gamma's in FCAL per event; counts/event", 14, -0.5, 13.5);
  center_axis(num_photons_cut);

  num_charged_cut = new TH1F("num_charged_cut", "; number of  charged tracks per event; counts/event", 14, -0.5, 13.5);
  center_axis(num_charged);
  //========================================================================//
  //>>>>>>>>>>>>> 2gamma histos            ///

  TDirectory * dir_2gamma = ROOTfile->mkdir("histos_2gamma");
  
  tree_2gamma = new TTree("tree_2gamma", "event number for background under 2gamma peak");
  tree_2gamma->SetDirectory(dir_2gamma);
  tree_2gamma->Branch("ev_num", &ev_num, "ev_num/I");

  mass_vs_elast_2gamma = new TH2F("mass_vs_elast_2gamma", "; E_{#eta} rec./E_{#gamma}; M_{#gamma#gamma} (GeV/c^{2})", 200, 0., 2., 1400, 0., 0.7);
  center_axis(mass_vs_elast_2gamma, dir_2gamma);

  mass_vs_total_e_2gamma = new TH2F("mass_vs_total_e_2gamma", "; E_{total, FCAL} (GeV); M_{#gamma#gamma} (GeV/c^{2})", 12000, 0., 12., 1400, 0., 0.7);
  center_axis(mass_vs_total_e_2gamma, dir_2gamma);

  theta_eta_2gamma = new TH1F("theta_eta_2gamma", "; #theta_{#eta, rec.} (deg.); counts/0.04#circ", 200, 0., 8.);
  center_axis(theta_eta_2gamma, dir_2gamma);

  mass_vs_elast_2gamma_cut = new TH2F("mass_vs_elast_2gamma_cut", "; E_{#eta} rec./E_{#gamma}; M_{#gamma#gamma} (GeV/c^{2})", 200, 0., 2., 1400, 0., 0.7);
  center_axis(mass_vs_elast_2gamma_cut, dir_2gamma);

  theta_eta_2gamma_cut = new TH1F("theta_eta_2gamma_cut", "; #theta_{#eta, rec.} (deg.); counts/0.04#circ", 200, 0., 8.);
  center_axis(theta_eta_2gamma_cut, dir_2gamma);

  //========================================================================//
  //>>>>>>>>>>>>>>> 6 gamma histos                    ///

  TDirectory * dir_6gamma = ROOTfile->mkdir("histos_6gamma");

  tree_6gamma = new TTree("tree_6gamma", "event number for background under 6gamma peak");
  tree_6gamma->SetDirectory(dir_6gamma);
  tree_6gamma->Branch("ev_num", &ev_num, "ev_num/I");

  mass_vs_elast_6gamma = new TH2F("mass_vs_elast_6gamma", "; E_{#eta} rec./E_{#gamma}; M_{#gamma#gamma} (GeV/c^{2})", 200, 0., 2., 1400, 0., 0.7);
  center_axis(mass_vs_elast_6gamma, dir_6gamma);

  mass_vs_total_e_6gamma = new TH2F("mass_vs_total_e_6gamma", "; E_{total, FCAL} (GeV); M_{#gamma#gamma} (GeV/c^{2})", 12000, 0., 12., 1400, 0., 0.7);
  center_axis(mass_vs_total_e_6gamma, dir_6gamma);

  theta_eta_6gamma = new TH1F("theta_eta_6gamma", "; #theta_{#eta, rec.} (deg.); counts/0.04#circ", 200, 0., 8.);
  center_axis(theta_eta_6gamma, dir_6gamma);

  mass_vs_elast_6gamma_cut = new TH2F("mass_vs_elast_6gamma_cut", "; E_{#eta} rec./E_{#gamma}; M_{#gamma#gamma} (GeV/c^{2})", 200, 0., 2., 1400, 0., 0.7);
  center_axis(mass_vs_elast_6gamma_cut, dir_6gamma);

  theta_eta_6gamma_cut = new TH1F("theta_eta_6gamma_cut", "; #theta_{#eta, rec.} (deg.); counts/0.04#circ", 200, 0., 8.);
  center_axis(theta_eta_6gamma_cut, dir_6gamma);
  
  //=========================================================================//
  fPDGdatabase = TDatabasePDG::Instance();

  return NOERROR;
}

//------------------------------------------------------------------
// evnt   -Fill histograms here
//------------------------------------------------------------------
jerror_t MyProcessor::evnt(JEventLoop *eventLoop, int eventnumber)
{
//   JEvent& event = eventLoop->GetJEvent();
//   JEventSource *source = event.GetJEventSource();
//   DEventSourceHDDM *hddm_source = dynamic_cast<DEventSourceHDDM*>(source);
//   if(!hddm_source){
//     cerr<<" This program MUST be used with an HDDM file as input!"<<endl;
//     exit(-1);
//   }
//   s_HDDM_t *hddm = (s_HDDM_t*)event.GetRef();
//   if(!hddm)return NOERROR;


//   cout<<"reaction type  = "<<(hddm->physicsEvents->in[0].reactions)->in[0].type<<endl;
//   cout<<"reaction weight  = "<<(hddm->physicsEvents->in[0].reactions)->in[0].weight<<endl;
//   cout<<"reaction mult  = "<<(hddm->physicsEvents->in[0].reactions)->mult<<endl;

//   cout<<"reaction run # = "<<hddm->physicsEvents->in[0].runNo<<endl;
//   cout<<"reaction evt # = "<<hddm->physicsEvents->in[0].eventNo<<endl;
//   cout<<"physics event multiplicity = "<<hddm->physicsEvents->mult<<endl<<endl;;

  // vector<const DPi0*> pions;
  vector<const DMCThrown*> throwns;
  vector<const DBeamPhoton*> beam_photon;
  // vector<const DFCALPhoton*> fcal_photons;
  vector<const DPhoton*> photons;
  //vector<const DParticle*> charged_part;

  eventLoop->Get(photons);
  // eventLoop->Get(pions);
  eventLoop->Get(throwns);
  eventLoop->Get(beam_photon);
  //eventLoop->Get(charged_part);

  pions.clear();
  status.clear();

  int N_FINAL_STATES =0;
  bool flag_proton  = false;
  bool flag_eta     = false;
  bool flag_charged = false;

  Reaction react;

  for(unsigned int i=0; i<throwns.size(); i++)
    {
      if(throwns[i]->mech==1)
	{
	  react.InsertFinalState(throwns[i]->pdgtype);
	  N_FINAL_STATES++;
	  if(throwns[i]->pdgtype==2212)flag_proton  = true;
	  else if(throwns[i]->pdgtype==221)flag_eta = true;
	  
	  if( throwns[i]->charge() !=0 && throwns[i]->pmag() > 0.14999 && throwns[i]->pdgtype!=2212)flag_charged = true;
	}
    }

  for(unsigned int i=0; i<beam_photon.size(); i++)
    E_beam->Fill(beam_photon[i]->lorentzMomentum().E());

  if(N_FINAL_STATES==2&&flag_eta==true&&flag_proton==true) return NOERROR;
  if(flag_charged == true) return NOERROR;




  if(photons.size()<2) return NOERROR;
  ev_num = eventnumber;
  double E_FCAL_Total = 0.;
  int N_FCAL_PHOTONS  = 0;
  int N_CHARGED_PART  = 0;  
  // bool DUMP_THROWNS = false;

  for(unsigned int i = 0; i < photons.size() ; i++)
    {
      if(photons[i]->getTag()==1)
	{
	  fcal_time->Fill(photons[i]->getTime()); 

	  DVector3 photon_pos = photons[i]->getPositionCal();
	  xy_position->Fill(photon_pos.X(), photon_pos.Y());

	  N_FCAL_PHOTONS++;
	  E_FCAL_Total += (photons[i]->lorentzMomentum()).E();
	}
      status.push_back(-1);
    }

  //N_CHARGED_PART = charged_part.size();

  num_photons->Fill(N_FCAL_PHOTONS);
  num_charged->Fill(N_CHARGED_PART);

  double E_beam_photon = beam_photon[0]->energy();

  total_energy_vs_elast->Fill(E_FCAL_Total/E_beam_photon, E_FCAL_Total);

    if(N_FCAL_PHOTONS<2) return NOERROR;
  //if(N_FCAL_PHOTONS!=2) return NOERROR;
  
  for (unsigned int i = 0; i < photons.size() ; i++) 
     {      
       //if(photons[i]->getTime()<-5 || photons[i]->getTime()>50)   continue;   
       // if (photons[i]->getdThetaCharge() < 0.05 ) continue;
       if(photons[i]->getTag() !=1 )continue;
       if(E_FCAL_Total>9.0)
	 {
	   fcal_time_cut->Fill(photons[i]->getTime());
	   num_photons_cut->Fill(N_FCAL_PHOTONS);
	   num_charged_cut->Fill(N_CHARGED_PART);
	   DVector3 photon_pos = photons[i]->getPositionCal();
	   xy_position_cut->Fill(photon_pos.X(), photon_pos.Y());
	 }
       for (unsigned int j = i+1; j < photons.size() ; j++) 
	 {
	   //if (photons[j]->getdThetaCharge() < 0.05 ) continue;
	   //  if(photons[j]->getTime()<-5 || photons[j]->getTime()>50)  continue;
	   
	   if(photons[j]->getTag() !=1) continue;
   
	   DLorentzVector gg_mom4 = photons[i]->lorentzMomentum() + photons[j]->lorentzMomentum();

	   double E_recon_2gamma = gg_mom4.E();
	   double ELAST_2gamma = E_recon_2gamma/E_beam_photon;
	   double theta_recon_2gamma = gg_mom4.Theta()*180./TMath::Pi();

	   double gg_mass = gg_mom4.M();

	   mass_vs_total_e_2gamma->Fill(E_FCAL_Total, gg_mass);
	   
	   if(E_FCAL_Total>9.0)  //in GeV
	     {
	       mass_vs_elast_2gamma->Fill(ELAST_2gamma, gg_mass);
	       theta_eta_2gamma->Fill(theta_recon_2gamma);
	       
	       if( gg_mass>0.48 && gg_mass<0.62 && ELAST_2gamma>0.89 && ELAST_2gamma<1.11)
		 {
		   //DUMP_THROWNS = true;
		   mass_vs_elast_2gamma_cut->Fill(ELAST_2gamma, gg_mass);
		   theta_eta_2gamma_cut->Fill(theta_recon_2gamma);
		   tree_2gamma->Fill();
		 } 
	     }

	   if(N_FCAL_PHOTONS>5 && status[i] ==-1 && status[j]==-1&& gg_mass>0.12 && gg_mass<0.165) 
	     {
	       pions.push_back(gg_mom4);
	       status[i]=1; status[j]=1;
	     }
	 }
     }

//   if(DUMP_THROWNS == true)
//     {
//       for(unsigned int i=0; i<throwns.size(); i++)
// 	{
// 	  cout<<"trackID = "  <<throwns[i]->myid<<" : "
// 	      <<"mechanism = "<<throwns[i]->mech<<" : "
// 	      <<"GEANT type = "<<throwns[i]->type<<" : "
// 	      <<"name = "<<fPDGdatabase->GetParticle(throwns[i]->pdgtype)->GetName()<<endl;
// 	}
//       cout<<endl;
//     }

  if(N_FCAL_PHOTONS<6 || pions.size()<3) return NOERROR;

  for(unsigned int i=0; i<pions.size(); i++)
    {
      for(unsigned int j=i+1; j<pions.size(); j++)
	{
	  for(unsigned int k=j+1; k<pions.size(); k++)
	    {
	      DLorentzVector mom4_6g = pions[i]+pions[j]+pions[k];
	      
	      double E_recon_6gamma = mom4_6g.E();
	      double ELAST_6gamma = E_recon_6gamma/E_beam_photon;
	      double theta_recon_6gamma = mom4_6g.Theta()*180./TMath::Pi();

	      double mass_6g = mom4_6g.M();

	      mass_vs_total_e_6gamma->Fill(mass_6g, E_FCAL_Total);

	      if(E_FCAL_Total>9.0) // in GeV
		{
		  mass_vs_elast_6gamma->Fill(ELAST_6gamma, mass_6g);
		  theta_eta_6gamma->Fill(theta_recon_6gamma);
	       
		  if( mass_6g>0.48 && mass_6g<0.62 && ELAST_6gamma>0.89 && ELAST_6gamma<1.11)
		    {
		      mass_vs_elast_6gamma_cut->Fill(ELAST_6gamma, mass_6g);
		      theta_eta_6gamma_cut->Fill(theta_recon_6gamma);
		      tree_6gamma->Fill();
		    } 
		}
	    }
	}
    }

//   if(eventnumber>140000) eventLoop->Quit();
  return NOERROR;
}

//------------------------------------------------------------------
// fini   -Close output file here
//------------------------------------------------------------------
jerror_t MyProcessor::fini(void)
{
	ROOTfile->Write();
	delete ROOTfile;
	cout<<endl<<"Closed ROOT file"<<endl;

	return NOERROR;
}