// $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 <PID/DPhoton.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"


//------------------
// 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 Photons");
	evt = new Event();
	evt_tree->Branch("E", &evt);

	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 DPhoton*> photons;
	vector<const DBCALPhoton*> bcalphotons;
	vector<const DFCALPhoton*> fcalphotons;

	loop->Get(photons);
	loop->Get(bcalphotons);
	loop->Get(fcalphotons);
	
	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);
	}

	// Event tree
	try{
		const DMCThrown *thrown;
		loop->GetSingle(thrown);
		evt->Clear();
		evt->thrown = thrown->lorentzMomentum();
		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 = 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();
		}
		evt_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;
}