// $Id$
//
//    File: DEventProcessor_splitoffs.cc
// Created: Tue Mar 16 16:05:43 EDT 2010
// Creator: davidl (on Darwin harriet.jlab.org 9.8.0 i386)
//

#include <TRACKING/DMCThrown.h>
#include <TRACKING/DMCTrajectoryPoint.h>
#include <PID/DPhoton.h>
#include <DLorentzVector.h>
#include <BCAL/DBCALHit.h>
#include <BCAL/DBCALTruthCell.h>
#include <BCAL/DBCALPhoton.h>
#include <FCAL/DFCALPhoton.h>

#include "DEventProcessor_splitoffs.h"
using namespace jana;


// Routine used to create our DEventProcessor
#include <JANA/JApplication.h>
extern "C"{
void InitPlugin(JApplication *app){
	InitJANAPlugin(app);
	app->AddProcessor(new DEventProcessor_splitoffs());
}
} // "C"

TLorentzVector MakeTlorentz(const DLorentzVector &src);


//------------------
// DEventProcessor_splitoffs (Constructor)
//------------------
DEventProcessor_splitoffs::DEventProcessor_splitoffs()
{

	pthread_mutex_init(&mutex, NULL);
}

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

}

//------------------
// init
//------------------
jerror_t DEventProcessor_splitoffs::init(void)
{
	// see Gluex-doc 1447
	dolbyc = new TH2D("dolbyc", "Energy Asymetry vs. opening angle", 100, 0.0, 2.0, 100, 0.0, 1.0);
	dolbyc->SetYTitle("1-A^{2}");
	dolbyc->SetXTitle("1-cos#psi");

	Nphoton = new TH1D("Nphoton","Number of reconstructed photons", 5, -0.5, 4.5);

	Nbcal_elements_per_cluster = new TH1D("Nbcal_elements_per_cluster", "Number of BCAL segments hit per cluster", 201, -0.5, 200.5);
	Nfcal_elements_per_cluster = (TH1D*)Nbcal_elements_per_cluster->Clone("Nfcal_elements_per_cluster");
	
	evt_tree = new TTree("event","Reconstructed Events");
	evt = new Event();
	evt_tree->Branch("A", &evt);

	reconPhoton_tree = new TTree("reconPhoton","Reconstructed Photons");
	reconPhoton = new ReconPhoton();
	reconPhoton_tree->Branch("R", &reconPhoton);

	bcalHit_tree = new TTree("bcalHit","BCAL Hits");
	bcalHit = new BCALHit();
	bcalHit_tree->Branch("H", &bcalHit);

	return NOERROR;
}

//------------------
// brun
//------------------
jerror_t DEventProcessor_splitoffs::brun(JEventLoop *eventLoop, int runnumber)
{
	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t DEventProcessor_splitoffs::evnt(JEventLoop *loop, int eventnumber)
{
	vector<const DBCALHit*> bcalhits;
	vector<const DBCALTruthCell*> truthcells;
	vector<const DPhoton*> photons;
	vector<const DBCALPhoton*> bcalphotons;
	vector<const DFCALPhoton*> fcalphotons;
	vector<const DMCTrajectoryPoint*> mctrajectorypoints;

	loop->Get(bcalhits);
	loop->Get(truthcells);
	loop->Get(photons);
	loop->Get(bcalphotons);
	loop->Get(fcalphotons);
	loop->Get(mctrajectorypoints);


	const DMCThrown *thrown;
	loop->GetSingle(thrown);
	float theta_gen = thrown->momentum().Theta();
	float z_gen = thrown->position().Z()+65.0/tan(theta_gen);
	float t_gen = 65.0/sin(theta_gen)/30.0;
	
	double sumTruth = 0.0;
	for(unsigned int i=0; i<truthcells.size(); i++)sumTruth += truthcells[i]->E;
	int Ntruthhits = (int)truthcells.size();
	
	// Find the last trajectory point corresponding to the primary particle
	const DMCTrajectoryPoint* traj_final=NULL;
	for(unsigned int i=0; i<mctrajectorypoints.size(); i++){
		if(mctrajectorypoints[i]->primary_track==1 && mctrajectorypoints[i]->track){
			traj_final = mctrajectorypoints[i];
		}
	}
	
	// Sort the BCAL hits into upstream and downstream lists. We'll need 
	// this when filling the BCALhits tree below, but this part can be done
	// outside the mutex lock so we do it here.
	vector<const DBCALHit*> uphits;
	vector<const DBCALHit*> downhits;
	for(unsigned int i=0; i<bcalhits.size(); i++){
		if(bcalhits[i]->end == DBCALGeometry::kUpstream)uphits.push_back(bcalhits[i]);
		else if(bcalhits[i]->end == DBCALGeometry::kDownstream)downhits.push_back(bcalhits[i]);		
	}

	// Find cells with both upstream and downstream hits and store
	// them in a pair. Again we do this for the BCALHit tree below,
	// but here since we are outside of the mutex lock.
	vector<pair<const DBCALHit*, const DBCALHit*> > bcal_coin;
	for(unsigned int i=0; i<uphits.size(); i++){
		const DBCALHit *uphit = uphits[i];
		for(unsigned int j=0; j<downhits.size(); j++){
			const DBCALHit *downhit = downhits[j];
			if(uphit->module != downhit->module)continue;
			if(uphit->layer  != downhit->layer )continue;
			if(uphit->sector != downhit->sector)continue;
			
			bcal_coin.push_back(pair<const DBCALHit*, const DBCALHit*>(uphit, downhit));
		}
	}
	
	// Lock mutex while we fill the tree and histos
	pthread_mutex_lock(&mutex);
	Nphoton->Fill(photons.size());
	
	// Dolby-C hist
	for(unsigned int i=0; i<photons.size(); i++){
		const DPhoton* photon1 = photons[i];
		for(unsigned int j=i+1; j<photons.size(); j++){
			const DPhoton* photon2 = photons[j];
			
			double E1 = photon1->energy();
			double E2 = photon2->energy();
			double A = (E1-E2)/(E1+E2);
			
			DVector3 mom1 = photon1->momentum();
			DVector3 mom2 = photon2->momentum();
			double psi = mom1.Angle(mom2);
			
			dolbyc->Fill(1-cos(psi), 1-A*A);
		}
	}
	
	// Loop over bcalphotons
	for(unsigned int i=0; i<bcalphotons.size(); i++){
		const DBCALPhoton *photon = bcalphotons[i];
		
		// The DBCALPhoton objects are made of one or more DBCALShower objects.
		// Each of those is made of one or more hits. 
		vector<const DBCALShower*> showers;
		photon->Get(showers);
		double Ncells_total = 0.0;
		for(unsigned int i=0; i<showers.size(); i++)Ncells_total+=(double)showers[i]->N_cell;
		
		Nbcal_elements_per_cluster->Fill(Ncells_total);
	}

	// Loop over fcalphotons
	for(unsigned int i=0; i<fcalphotons.size(); i++){
		const DFCALPhoton *photon = fcalphotons[i];
		
		// The DFCALPhoton objects are made of one or more DBCALShower objects.
		// Each of those is made of one or more hits. 
		vector<const DFCALCluster*> clusters;
		photon->Get(clusters);
		double Ncells_total = 0.0;
		for(unsigned int i=0; i<clusters.size(); i++)Ncells_total+=(double)clusters[i]->getHits();
		
		Nfcal_elements_per_cluster->Fill(Ncells_total);
	}

	try{
		// Event tree
		evt->Clear();
		evt->event = eventnumber;
		evt->thrown = MakeTlorentz(thrown->lorentzMomentum());
		evt->sumTruth = sumTruth;
		evt->Ntruthhits = Ntruthhits;
		for(unsigned int i=0; i<photons.size(); i++){
			const DPhoton* photon = photons[i];
			TClonesArray &prts = *(evt->recon);
			Particle *prt = new(prts[evt->Nrecon++]) Particle();
			prt->p = MakeTlorentz(photon->lorentzMomentum());
			double R = 65.0;
			double drphi = R*evt->thrown.DeltaPhi(prt->p);
			double dz = evt->thrown.Z() - prt->p.Z();
			prt->dist = sqrt(drphi*drphi + dz*dz);
			prt->tag = photon->getTag();
			prt->idx = -1; // fill in after sorting

			// Get Shower object used to create this DPhoton so we can get
			// the value of E_raw. (If either GetSingle() calls fail below
			// an exception should be thrown).
			try{
				const DBCALPhoton *bcalphoton;
				const DBCALShower *bcalshower;
				photon->GetSingle(bcalphoton);
				bcalphoton->GetSingle(bcalshower);
				prt->E_raw = bcalshower->E_raw;
			}catch(...){}
		}
		
		// Sort by energy decreasing
		evt->recon->Sort();
		for(UInt_t i=0; i<evt->Nrecon; i++)((Particle*)evt->recon->At(i))->idx = i;

		if(traj_final){
			evt->x = traj_final->x;
			evt->y = traj_final->y;
			evt->z = traj_final->z;
			evt->t = traj_final->t;
			evt->px = traj_final->px;
			evt->py = traj_final->py;
			evt->pz = traj_final->pz;
			evt->E = traj_final->E;
			evt->part = traj_final->part;
			evt->mech = traj_final->mech;
		}
		
		evt_tree->Fill();
		
		// ReconPhoton tree

		// First, we have to make a pass over all reconstructed photons in
		// order to make a list of the distance each is from the thrown projection.
		// This is so we can record the position of each reconstructed photon
		// in that list when it is ordered by distance.
		multimap<double,int> dist_to_throwns; // key=distance value=index
		for(unsigned int i=0; i<photons.size(); i++){
			const DPhoton* photon = photons[i];
			DVector3 dpos_cal = photon->getPositionCal();
			TVector3 pos_cal(dpos_cal.X(), dpos_cal.Y(), dpos_cal.Z());

			// thrown particle's projection on BCAL surface
			TVector3 pos_thrown = evt->thrown.Vect();
			pos_thrown *= 65.0/pos_thrown.Perp();
			pos_thrown += TVector3(0.0, 0.0, 65.0); // shift to center of target

			double dist_to_thrown = (pos_cal-pos_thrown).Mag();
			dist_to_throwns.insert(pair<double,int>(dist_to_thrown, i));
		}
		
		// Now, we loop over reconstructed photons a second time, creating
		// an entry in the reconPhoton tree for each
		for(unsigned int i=0; i<photons.size(); i++){
			const DPhoton* photon = photons[i];
			DVector3 dpos_cal = photon->getPositionCal();
			
			// For the moment, we are only interested in BCAL
			if(dpos_cal.Z() > 600)continue;
			
			// Find the distance to the thrown particle's projection on surface
			// of calorimeter. This a bit complicated due to the idist_to_thrown
			// field needing to know the position of this in a list ordered
			// by dist_to_thrown. Those are all calculated in the loop above. Here,
			// we need to identify the position of the current reconstructed
			// photon in the list.
			multimap<double,int>::iterator iter;
			int idist_to_thrown=0;
			double dist_to_thrown = 0.0;
			for(iter=dist_to_throwns.begin(); iter!=dist_to_throwns.end(); iter++, idist_to_thrown++){
				if(iter->second == (int)i){
					dist_to_thrown = iter->first;
					break;
				}
			}
			
			reconPhoton->Clear();
			reconPhoton->event = eventnumber;
			reconPhoton->thrown = evt->thrown;
			reconPhoton->recon = MakeTlorentz(photon->lorentzMomentum());
			reconPhoton->pos.SetXYZ(dpos_cal.X(), dpos_cal.Y(), dpos_cal.Z());
			reconPhoton->t = photon->getTime();
			reconPhoton->Nrecon = (int)photons.size();
			reconPhoton->irecon = (int)i;
			reconPhoton->dist_to_thrown = dist_to_thrown;
			reconPhoton->idist_to_thrown = idist_to_thrown;
			reconPhoton->theta_gen = theta_gen;
			reconPhoton->z_gen     = z_gen;
			reconPhoton->t_gen     = t_gen;
			
			const DMCTrajectoryPoint *traj = GetMaxETrajectoryPoint(mctrajectorypoints, reconPhoton->pos, 15.0/* cm */);
			if(traj){
				reconPhoton->ptype_closest = traj->part;
				reconPhoton->pos_closest = TVector3(traj->x, traj->y, traj->z);
				traj = GetFirstTrajectoryPoint(mctrajectorypoints, traj);
				if(traj){
					reconPhoton->pos_origin = TVector3(traj->x, traj->y, traj->z);
					reconPhoton->mech = traj->mech;
				}
			}

			reconPhoton_tree->Fill();
		}
		
		// BCALHit Tree
		// The BCAL hit tree keeps a single entry for every cell that has
		// both ends hit. Assume each cell has at most one hit per end.
		// "Coincidences" (i.e. both ends of cell hit) are found above
		// before the mutex gets locked. Loop over them here and use them
		// to fill  the BCALHit tree.
		for(unsigned int i=0; i<bcal_coin.size(); i++){
			const DBCALHit *uphit = bcal_coin[i].first;
			const DBCALHit *downhit = bcal_coin[i].second;
			
			bcalHit->Clear();
			
			bcalHit->event     = eventnumber;
			bcalHit->tup       = uphit->t;
			bcalHit->tdown     = downhit->t;
			bcalHit->Eup       = uphit->E;
			bcalHit->Edown     = downhit->E;
			bcalHit->module    = uphit->module;
			bcalHit->layer     = uphit->layer;
			bcalHit->sector    = uphit->sector;
			bcalHit->theta_gen = theta_gen;
			bcalHit->z_gen     = z_gen;
			bcalHit->t_gen     = t_gen;

			bcalHit_tree->Fill();
		}
		
		
	}catch(...){}

	pthread_mutex_unlock(&mutex);

	return NOERROR;
}

//------------------
// erun
//------------------
jerror_t DEventProcessor_splitoffs::erun(void)
{
	// Any final calculations on histograms (like dividing them)
	// should be done here. This may get called more than once.
	return NOERROR;
}

//------------------
// fini
//------------------
jerror_t DEventProcessor_splitoffs::fini(void)
{
	// Called at very end. This will be called only once

	return NOERROR;
}

//------------------
// GetMaxETrajectoryPoint
//------------------
const DMCTrajectoryPoint* DEventProcessor_splitoffs::GetMaxETrajectoryPoint(vector<const DMCTrajectoryPoint*> mctrajectorypoints, TVector3 &pos, float max_dist)
{
	// Loop over all trajectory points and find the one that has the most 
	// energy that is within mas_dist of pos.

	const DMCTrajectoryPoint *traj_max = NULL;
	double Emax = 0.0;
	for(unsigned int i=0; i<mctrajectorypoints.size(); i++){
		const DMCTrajectoryPoint *traj = mctrajectorypoints[i];
		
		// GEANT will create some massive secondary particles by simply
		// adding kinetic energy and assuming the particle was
		// already there as part of the material. It does this with
		// protons and neutrons, but also with Deuterons, Tritium,
		// and Alpha particles.For these, the available energy is
		// just the kinetic energy of the particle while for things 
		// like pions, the total energy should be included. We handle
		// this by setting a value for the "inaccessible" energy that
		// is tied up in the mass. For pions, we consider this zero.
		double unavailable_mass = 0.0;
		switch(traj->part){
			case 13:	unavailable_mass= 0.93956563;   break;
			case 14:	unavailable_mass= 0.93827231;   break;
			case 45:	unavailable_mass= 1.87561300;   break;
			case 46:	unavailable_mass= 2.80925000;   break;
			case 47:	unavailable_mass= 3.72741700;   break;
		}
		double Eavailable = traj->E - unavailable_mass;
		
		TVector3 pos_traj(traj->x, traj->y, traj->z);
		if((pos_traj-pos).Mag() <= max_dist){
			if(traj_max != NULL){
				if(Eavailable > Emax){
					traj_max = traj;
					Emax = Eavailable;
				}
			}else{
				traj_max = traj;
			}
		}
	}
	
	return traj_max;
}

//------------------
// GetFirstTrajectoryPoint
//------------------
const DMCTrajectoryPoint* DEventProcessor_splitoffs::GetFirstTrajectoryPoint(vector<const DMCTrajectoryPoint*> mctrajectorypoints, const DMCTrajectoryPoint *traj)
{
	// Loop over all trajectory points and find the first in the list
	// that matches the given trajectory point. We assume the order is
	// maintained from what GEANT produced making the first in the list
	// the birth point.
	const DMCTrajectoryPoint *traj_origin = traj;
	for(unsigned int i=0; i<mctrajectorypoints.size(); i++){
		const DMCTrajectoryPoint *my_traj = mctrajectorypoints[i];
		
		if(my_traj->track != traj->track)continue;
		if(my_traj->primary_track != traj->primary_track)continue; // redundant
		if(my_traj->part != traj->part)continue;	// doubly redundant
		
		traj_origin = my_traj;
		break;
	}

	return traj_origin;
}


//------------------
// MakeTlorentz
//------------------
TLorentzVector MakeTlorentz(const DLorentzVector &src)
{
	return TLorentzVector(src.X(), src.Y(), src.Z(), src.E());
}