// $Id$
//
//    File: JEventProcessor_ebode_fcal.cc
// Created: Fri  7 Jul 11:49:30 EDT 2017
// Creator: ebode (on Linux ifarm1402.jlab.org 3.10.0-327.el7.x86_64 x86_64)
//

#include "JEventProcessor_ebode_fcal.h"
using namespace jana;

#include <particleType.h>
#include <DANA/DApplication.h>
#include <HDGEOMETRY/DGeometry.h>
#include <FCAL/DFCALShower.h>
#include <TRACKING/DTrackTimeBased.h>
#include <FCAL/DFCALHit.h>
#include <FCAL/DFCALGeometry.h>

// Routine used to create our JEventProcessor
#include <JANA/JApplication.h>
#include <JANA/JFactory.h>

extern "C"{
void InitPlugin(JApplication *app){
	InitJANAPlugin(app);
	app->AddProcessor(new JEventProcessor_ebode_fcal());
}
} // "C"


//------------------
// JEventProcessor_ebode_fcal (Constructor)
//------------------
JEventProcessor_ebode_fcal::JEventProcessor_ebode_fcal()
{
	zFCAL_front = 0.0;
	zFCAL_back = 0.0;
}

//------------------
// ~JEventProcessor_ebode_fcal (Destructor)
//------------------
JEventProcessor_ebode_fcal::~JEventProcessor_ebode_fcal()
{

}

//------------------
// init
//------------------
jerror_t JEventProcessor_ebode_fcal::init(void)
{
	hfcalShowerE = new TH1D("hfcalShowerE", "FCAL shower Energy;Energy (GeV)", 200, 0.0, 8.0); 
	hfcalShowerP = new TH2D("hfcalShowerP", "FCAL shower Position;Position", 280, -140.0, 140.0, 280, -140.0, 140.0);
	hfcalShowerT = new TH1D("hfcalShowerT", "FCAL shower Time;Time", 120, -10.0, 50.0);	
	
	tPionTrack = new TTree("tPionTrack", "Charged pion tracks hitting FCAL");
	tPionTrack->Branch("t", "pion_track", &pi_trk);
	
	return NOERROR;
}

//------------------
// brun
//------------------
jerror_t JEventProcessor_ebode_fcal::brun(JEventLoop *eventLoop, int32_t runnumber)
{
	DApplication *dapp = dynamic_cast<DApplication*>( eventLoop->GetJApplication() );
	const DGeometry *dgeom = dapp->GetDGeometry(runnumber);
	dgeom->GetFCALZ(zFCAL_front);
	zFCAL_back = zFCAL_front + 45.0;
	jout << " FCAL z-front: " << zFCAL_front << endl;
	jout << " FCAL  z-back: " << zFCAL_back << endl;

	// This is called whenever the run number changes
	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t JEventProcessor_ebode_fcal::evnt(JEventLoop *loop, uint64_t eventnumber)
{
	vector<const DFCALShower*> fcalshowers;
	vector<const DTrackTimeBased*> tbts;
	vector<const DFCALHit*> fcalhit;
	vector<const DFCALGeometry*> fcalgeometries;
	loop->Get(fcalshowers);
	loop->Get(tbts);
	loop->Get(fcalhit);
	loop->Get(fcalgeometries);
	
	// Get DFCALGeometry object
	const DFCALGeometry* fcalgeom = fcalgeometries.empty() ? NULL:fcalgeometries[0];
	if( fcalgeom==NULL ){
		jerr << "NO DFCALGeometry OBJECTS!" << endl;
		exit(-1);
	}
	
	map<const DTrackTimeBased*, const DFCALShower*> trk_fcalshower;	
	map<const DTrackTimeBased*, const DFCALHit*   > trk_fcalhit;
	map<const DTrackTimeBased*, DVector3> pos_front;
	map<const DTrackTimeBased*, DVector3> pos_back;
	map<const DTrackTimeBased*, vector<const DFCALHit*> > hits9;
	map<const DTrackTimeBased*, vector<const DFCALHit*> > hits25;
	for(auto tbt : tbts){
		auto rt = tbt->rt;
		
		if( (tbt->PID()!=PiPlus) && (tbt->PID()!=PiMinus) ) continue;
		
		// Project track to front face of FCAL
		DVector3 origin(0.0, 0.0, zFCAL_front);
		DVector3 norm(0.0, 0.0, -1.0);
		DVector3 pos;
		if( NOERROR == rt->GetIntersectionWithPlane(origin, norm, pos)){
			pos_front[tbt] = pos;
			
			// Search for blocks closest to track projection
			int row = fcalgeom->row((float)pos.Y());
			int col = fcalgeom->column((float)pos.X());
			for(auto hit : fcalhit ){
				if( abs(hit->row - row)<=1 && abs(hit->column - col)<=1 ) hits9[tbt].push_back(hit);
				if( abs(hit->row - row)<=2 && abs(hit->column - col)<=2 ) hits25[tbt].push_back(hit);
			}

			// Match tracks to FCAL hit and shower
			for(auto shower : fcalshowers ){
				double r = (pos - shower->getPosition()).Perp();
				if( r<5.0 ) trk_fcalshower[tbt] = shower;
			}
			double Emax = 0.0;
			for(auto hit : hits25[tbt]){
				if(hit->E > Emax){
					Emax = hit->E;
					trk_fcalhit[tbt] = hit;
				}
			}
		}

		// Project track to back face of FCAL
		origin.SetXYZ(0.0, 0.0, zFCAL_back);
		if( NOERROR == rt->GetIntersectionWithPlane(origin, norm, pos)){			
			pos_back[tbt] = pos;
		}
		
	}
	
	// -------------------------------------------------------------------------
	japp->RootFillLock(this);
	
	// Fill pion track tree
	for(auto tbt : tbts){
	
		if( (tbt->PID()!=PiPlus) && (tbt->PID()!=PiMinus) ) continue;
		
		double E9  = 0.0;
		double E25 = 0.0;
		for(auto hit : hits9[tbt] ) E9 += hit->E;
		for(auto hit : hits25[tbt]) E25 += hit->E;
		
		double r_front = -1.0;
		double r_back  = -1.0;
		if(pos_front.count(tbt)+trk_fcalhit.count(tbt) == 2){
			double dx = pos_front[tbt].X() - trk_fcalhit[tbt]->x;
			double dy = pos_front[tbt].Y() - trk_fcalhit[tbt]->y;
			r_front = sqrt( dx*dx + dy*dy );
		}
		if(pos_back.count(tbt)+trk_fcalhit.count(tbt) == 2){
			double dx = pos_back[tbt].X() - trk_fcalhit[tbt]->x;
			double dy = pos_back[tbt].Y() - trk_fcalhit[tbt]->y;
			r_back = sqrt( dx*dx + dy*dy );
		}

		pi_trk.ptype = tbt->PID();
		pi_trk.Etrk = tbt->energy();
		pi_trk.Ehit = (trk_fcalhit[tbt]==NULL ? 0.0:trk_fcalhit[tbt]->E);
		pi_trk.E9  = E9;
		pi_trk.E25 = E25;
		pi_trk.Eshower = (trk_fcalshower[tbt]==NULL ? 0.0:trk_fcalshower[tbt]->getEnergy());
		pi_trk.N9 = hits9[tbt].size();
		pi_trk.N25 = hits25[tbt].size();
		pi_trk.p = tbt->momentum().Mag();
		pi_trk.theta = tbt->momentum().Theta();
		pi_trk.phi = tbt->momentum().Phi();
		pi_trk.r_front = r_front;
		pi_trk.r_back = r_back;
		
		if(pos_front.count(tbt)==1){
			pi_trk.pos_x_front = pos_front[tbt].X();
			pi_trk.pos_y_front = pos_front[tbt].Y();
		}
		if(pos_back.count(tbt)==1){
			pi_trk.pos_x_back = pos_back[tbt].X();
			pi_trk.pos_y_back = pos_back[tbt].Y();
		}
		if(trk_fcalhit.count(tbt)==1 && trk_fcalhit[tbt]!=NULL){
			pi_trk.pos_x_block = trk_fcalhit[tbt]->x;
			pi_trk.pos_y_block = trk_fcalhit[tbt]->y;
		}
		
		tPionTrack->Fill();
	}	


	// FCAL Shower histograms
	for(auto shower : fcalshowers ){
		hfcalShowerE->Fill( shower->getEnergy()); 
		hfcalShowerP->Fill( shower->getPosition().x(),shower->getPosition().y());
		hfcalShowerT->Fill( shower->getTime());
	}


	japp->RootFillUnLock(this); //RELEASE ROOT LOCK
	// -------------------------------------------------------------------------


	return NOERROR;
}

//------------------
// erun
//------------------
jerror_t JEventProcessor_ebode_fcal::erun(void)
{
	// This is called whenever the run number changes, before it is
	// changed to give you a chance to clean up before processing
	// events from the next run number.
	return NOERROR;
}

//------------------
// fini
//------------------
jerror_t JEventProcessor_ebode_fcal::fini(void)
{
	// Called before program exit after event processing is finished.
	return NOERROR;
}