// $Id$
//
//    File: JEventProcessor_fcaldave.cc
// Created: Sat Nov 22 22:47:43 EST 2014
// Creator: davidl (on Linux gluon47.jlab.org 2.6.32-358.23.2.el6.x86_64 x86_64)
//

#include <iostream>
using namespace std;

#include "JEventProcessor_fcaldave.h"
using namespace jana;


#include "DFCALDigiHit_factory_DAVECUT.h"
#include "JFactoryGenerator_DFCALDigiHit_DAVECUT.h"

//double E_CALIB = 0.135/0.038;
double E_CALIB = 1.0;

// 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_fcaldave());
	app->AddFactoryGenerator(new JFactoryGenerator_DFCALDigiHit_DAVECUT());
}
} // "C"


//------------------
// JEventProcessor_fcaldave (Constructor)
//------------------
JEventProcessor_fcaldave::JEventProcessor_fcaldave()
{


}

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

}

//------------------
// init
//------------------
jerror_t JEventProcessor_fcaldave::init(void)
{
	japp->RootWriteLock();

	ped_vs_id = new TH2D("ped_vs_id", "FCAL pedestal vs. id", 2800, 0.5, 2800.5, 200, 0.0, 200.0);
	amp_vs_id = new TH2D("amp_vs_id", "FCAL Amplitude vs. id", 2800, 0.5, 2800.5, 200, -10.0, 20000.0);
	tdc_vs_id = new TH2D("tdc_vs_id", "FCAL TDC vs. id", 2800, 0.5, 2800.5, 200, 0.0, 4000.0);
	integral_vs_peak = new TH2D("integral_vs_peak", "FCAL integral vs. peak", 200, 0.0, 4000.0, 200, 0.0, 40000.0);
	Ncluster = new TH1D("Ncluster", "Number of clusters", 21, -0.5, 20.5);

	E_vs_id   = new TH2D("E_vs_id", "FCAL Energy vs. id", 2800, 0.5, 2800.5, 205, -0.5, 10.0);
	t_vs_id = new TH2D("t_vs_id", "FCAL time vs. id", 2800, 0.5, 2800.5, 400, 0.0, 400.0);
	E_cluster = new TH1D("E_cluster", "FCAL Cluster energy", 1050, -0.5, 10.0);
	t_cluster = new TH1D("t_cluster", "FCAL Cluster time", 1100, -100.0, 1000.0);
	E_cluster2 = new TH1D("E_cluster2", "FCAL Cluster energy when 2 clusters", 1050, -0.5, 10.0);
	E_cluster_time_cut = new TH1D("E_cluster_time_cut", "FCAL Cluster energy cut on time 86<=t<=133", 1050, -0.5, 10.0);
	xy_cluster = new TH2D("xy_cluster", "FCAL cluster x/y", 240, -120.0, 120.0, 240, -120.0, 120.0);
	z_cluster = new TH1D("z_cluster", "FCAL Cluster z", 1000, 0.0, 1000.0);


	mass = new TH1D("mass", "FCAL Cluster 2 clusters invariant mass", 2500, -0.5, 2.0);
	z_pi0 = new TH1D("z_pi0", "z of pi0 decay vertex assuming pi0 mass", 2400, -200.0, 1000.0);

	pi0_ptr = & pi0;
	pi0_tree = new TTree("pi0", "Two photon cluster combos");
	pi0_tree->Branch("T","pi0_t", &pi0_ptr);

	japp->RootUnLock();

	return NOERROR;
}

//------------------
// brun
//------------------
jerror_t JEventProcessor_fcaldave::brun(JEventLoop *eventLoop, int runnumber)
{
	// This is called whenever the run number changes
	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t JEventProcessor_fcaldave::evnt(JEventLoop *loop, int eventnumber)
{
	vector<const DFCALDigiHit*> fcaldigihits;
	vector<const DFCALHit*> fcalhits;
	vector<const DFCALCluster*> fcalclusters;
	vector<const DTrackWireBased*> trackwirebased;
	vector<const DVertex*> vertexes;
	loop->Get(fcaldigihits);
	loop->Get(fcalhits);
	loop->Get(fcalclusters);
	loop->Get(trackwirebased);
	loop->Get(vertexes);

	// Get track projections onto FCAL
	vector<track_projection> tprojs;
	GetProjections(trackwirebased, tprojs);

	japp->RootWriteLock();
	
	// DFCALDigiHit
	for(uint32_t i=0; i<fcaldigihits.size(); i++){
		const DFCALDigiHit *fcaldigihit = fcaldigihits[i];
		
		uint32_t pedestal = 101;
		const Df250PulsePedestal *pp = NULL;
		fcaldigihit->GetSingle(pp);
		if(pp) pedestal = pp->pedestal;
		if(pp) integral_vs_peak->Fill((double)pp->pulse_peak, (double)fcaldigihit->pulse_integral);
		
		uint32_t id = fcalgeom.channel(fcaldigihit->row, fcaldigihit->column);
		ped_vs_id->Fill(id, pedestal);
		tdc_vs_id->Fill(id, fcaldigihit->pulse_time);
		
		if(pedestal<97) pedestal = 101;
		
		uint32_t Nsamples = fcaldigihit->nsamples_integral;
		
		double A = (double)fcaldigihit->pulse_integral - (double)(pedestal*Nsamples);
		amp_vs_id->Fill(id, A);
		
	}
	
	// DFCALHit
	double Efcal = 0.0;
	for(uint32_t i=0; i<fcalhits.size(); i++){
		const DFCALHit *fcalhit = fcalhits[i];

		uint32_t id = fcalgeom.channel(fcalhit->row, fcalhit->column);

		E_vs_id->Fill(id, fcalhit->E);
		t_vs_id->Fill(id, fcalhit->t);
		
		Efcal += fcalhit->E;
	}

	// DFCALCluster
	Ncluster->Fill((double)fcalclusters.size());
	for(uint32_t i=0; i<fcalclusters.size(); i++){
		const DFCALCluster *fcalcluster = fcalclusters[i];
		
		double E = fcalcluster->getEnergy();
		double t = fcalcluster->getTime();

		E_cluster->Fill(E);
		t_cluster->Fill(t);
		if(fcalclusters.size()>1) E_cluster2->Fill(E);
		if(t>=86 && t<=133) E_cluster_time_cut->Fill(E);
		
		TVector3 pos = fcalcluster->getCentroid();
		xy_cluster->Fill(pos.X(), pos.Y());
		z_cluster->Fill(pos.Z());
		
	}
	if(fcalclusters.size()>1){
		for(uint32_t i=0; i<fcalclusters.size()-1; i++){
			const DFCALCluster *fcalcluster1 = fcalclusters[i];
			double E1 = fcalcluster1->getEnergy();
			E1 *= E_CALIB;
			TVector3 pos1 = fcalcluster1->getCentroid();
			TVector3 mom1 = pos1;
			mom1.SetMag(E1);
			TLorentzVector p1(mom1, E1);

			float dist_track1;
			track_projection *proj1 = FindClosestTrack(tprojs, fcalcluster1, dist_track1);

			for(uint32_t j=i+1; j<fcalclusters.size(); j++){
				const DFCALCluster *fcalcluster2 = fcalclusters[j];
				double E2 = fcalcluster2->getEnergy();
				E2 *= E_CALIB;
				TVector3 pos2 = fcalcluster2->getCentroid();
				TVector3 mom2 = pos2;
				mom2.SetMag(E2);
				TLorentzVector p2(mom2, E2);

				float dist_track2;
				track_projection *proj2 = FindClosestTrack(tprojs, fcalcluster2, dist_track2);

			if(E1>0.2 && E2>0.2){
					double m = (p1+p2).M();
					if(m > 0.010) mass->Fill(m);
					
					// Calculate z of the pi0 vertex assuming a pi0 mass
					double mpi = 0.135;
					double C = 1.0 - mpi*mpi/(2.0*E1*E2);
					double C2 = C*C;
					double r1sq = pos1.Perp2();
					double r2sq = pos2.Perp2();
					double alpha = pos1.X()*pos2.X() + pos1.Y()*pos2.Y();
					double a = 1.0 - C2;
					double b = 2.0*alpha - (r1sq+r2sq)*C2;
					double c = alpha*alpha - C2*r1sq*r2sq;
					
					double disc = sqrt(b*b - 4.0*a*c);
					double z2  = (-b + disc)/(2.0*a);
					double z = isfinite(z2) ? sqrt( z2 ):0.0;
					if( z != 0.0)z_pi0->Fill(z);
					

					
					// Fill tree
					pi0.Nhits1 = fcalcluster1->getHits();
					pi0.Nhits2 = fcalcluster2->getHits();
					pi0.Nclusters = (int)fcalclusters.size();
					pi0.Nwbtracks = (int)trackwirebased.size();
					pi0.E1 = E1;
					pi0.E2 = E2;
					pi0.r1 = sqrt(r1sq);
					pi0.r2 = sqrt(r2sq);
					pi0.t1 = fcalcluster1->getTime();
					pi0.t2 = fcalcluster2->getTime();
					pi0.dist_track1 = dist_track1;
					pi0.dist_track2 = dist_track2;
					pi0.p_track1 = proj1!=NULL ? proj1->p:0.0;
					pi0.p_track2 = proj2!=NULL ? proj2->p:0.0;
					pi0.Efcal = Efcal;
					pi0.z_pi0 = z;
					pi0.mass = m;
					pi0.mom1 = mom1;
					pi0.mom2 = mom2;
					pi0.pos1 = pos1;
					pi0.pos2 = pos2;
					pi0.vertex = vertexes.empty() ? TVector3(1000.0, 1000.0, 1000.0):vertexes[0]->dSpacetimeVertex.Vect();
					pi0_tree->Fill();
				}
			}
		}
	}

	japp->RootUnLock();


	return NOERROR;
}

//------------------
// GetProjections
//------------------
uint32_t JEventProcessor_fcaldave::GetProjections(vector<const DTrackWireBased*> &wbtrks, vector<track_projection> &tprojs)
{
	// Location of front face of FCAL
	DVector3 fcalpos(0.0, 0.0, 560.275 + 65.0);
	
	for(uint32_t i=0; i<wbtrks.size(); i++){
		const DTrackWireBased *trk = wbtrks[i];
		
		DVector3 norm(0.0,0.0,-1);
		DVector3 pos,mom;
		if(trk->rt->GetIntersectionWithPlane(fcalpos,norm,pos,mom)==NOERROR){

			track_projection tproj;
			tproj.x = pos.X();
			tproj.y = pos.Y();
			tproj.p = mom.Mag();
			tproj.trk = trk;
		
			tprojs.push_back(tproj);
		}
	}
	
	return tprojs.size();
}

//------------------
// FindClosestTrack
//------------------
JEventProcessor_fcaldave::track_projection* JEventProcessor_fcaldave::FindClosestTrack(vector<track_projection> &tprojs, const DFCALCluster *cluster, float &dist)
{
	TVector3 pos = cluster->getCentroid();
	float x_cluster = pos.X();
	float y_cluster = pos.Y();

	track_projection *tproj_best = NULL;
	float d2_min = 1.0E6;
	for(uint32_t i=0; i<tprojs.size(); i++){
		track_projection *tproj = &tprojs[i];
		
		float dx = x_cluster - tproj->x;
		float dy = y_cluster - tproj->y;
		float d2 = dx*dx - dy*dy;
		
		if(d2<d2_min || tproj_best==NULL){
			tproj_best = tproj;
			d2_min = d2;
		}
	}
	
	dist = sqrt(d2_min);
	
	return tproj_best;
}

//------------------
// erun
//------------------
jerror_t JEventProcessor_fcaldave::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_fcaldave::fini(void)
{
	// Called before program exit after event processing is finished.
	return NOERROR;
}