// $Id$
//
//    File: DEventProcessor_bcal_calib.cc
// Created: Fri Jun  3 00:15:38 EDT 2011
// Creator: davidl (on Linux ifarml1 2.6.18-128.el5)
//

#include <TRACKING/DMCThrown.h>
#include <DLorentzVector.h>
#include <BCAL/DBCALHit.h>
#include <BCAL/DBCALShower.h>

#include "DEventProcessor_bcal_calib.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_bcal_calib());
}
} // "C"

TLorentzVector MakeTlorentz(const DLorentzVector &src);


//------------------
// DEventProcessor_bcal_calib (Constructor)
//------------------
DEventProcessor_bcal_calib::DEventProcessor_bcal_calib()
{

	pthread_mutex_init(&mutex, NULL);
}

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

}

//------------------
// init
//------------------
jerror_t DEventProcessor_bcal_calib::init(void)
{

	Nshower = new TH1D("Nshower","Number of reconstructed showers in BCAL", 21, -0.5, 20.5);
	Nhits_per_shower = new TH1D("Nhits_per_shower", "Number of BCAL hits per shower", 201, -0.5, 200.5);
	
	reconPhoton_tree = new TTree("reconPhoton","Reconstructed Photons");
	reconPhoton = new ReconPhoton();
	reconPhoton_tree->Branch("R", &reconPhoton);

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

	return NOERROR;
}

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

//------------------
// evnt
//------------------
jerror_t DEventProcessor_bcal_calib::evnt(JEventLoop *loop, int eventnumber)
{
	const DMCThrown *thrown;
	vector<const DBCALHit*> bcalhits;
	vector<const DBCALShower*> bcalshowers;

	loop->GetSingle(thrown);
	loop->Get(bcalhits);
	loop->Get(bcalshowers, "KLOE");
	
	// 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);
	Nshower->Fill(bcalshowers.size());
	
	try{		

		// ReconPhoton tree

		// Loop over bcalshowers
		for(unsigned int i=0; i<bcalshowers.size(); i++){
			const DBCALShower *shower = bcalshowers[i];
		
			Nhits_per_shower->Fill((double)shower->N_cell);
		
			reconPhoton->Clear();
			reconPhoton->event = eventnumber;
			reconPhoton->E = shower->E;
			reconPhoton->E_raw = shower->E_raw;
			reconPhoton->pos.SetXYZ(shower->x, shower->y, shower->z);
			reconPhoton->t = shower->t;
			reconPhoton->N_cell = shower->N_cell;
			reconPhoton->Nrecon = (int)bcalshowers.size();
			reconPhoton->E_thrown = thrown->energy();
			reconPhoton->p_thrown = thrown->momentum();
		
			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 = thrown->momentum().Theta();
			bcalHit->z_gen     = thrown->position().Z();
			bcalHit->t_gen     = thrown->t0();

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

	pthread_mutex_unlock(&mutex);

	return NOERROR;
}

//------------------
// erun
//------------------
jerror_t DEventProcessor_bcal_calib::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_bcal_calib::fini(void)
{
	// Called at very end. This will be called only once

	return NOERROR;
}

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