// $Id$
//
//    File: JEventProcessor_bcal_leakage.cc
// Created: Mon Oct  3 05:36:10 EDT 2011
// Creator: davidl (on Linux ifarm1101 2.6.18-128.7.1.el5 x86_64)
//

#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;

#include <TRACKING/DMCTrajectoryPoint.h>
#include <TRACKING/DMCThrown.h>

#include "JEventProcessor_bcal_leakage.h"
using namespace jana;


ClassImp(leakhit);

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


//------------------
// JEventProcessor_bcal_leakage (Constructor)
//------------------
JEventProcessor_bcal_leakage::JEventProcessor_bcal_leakage()
{

}

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

}

//------------------
// init
//------------------
jerror_t JEventProcessor_bcal_leakage::init(void)
{
	lhitptr = &lhit;

	leaks = new TTree("leaks", "BCAL energy leakage");
	leaks->Branch("L",&lhit);

	return NOERROR;
}

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

//------------------
// evnt
//------------------
jerror_t JEventProcessor_bcal_leakage::evnt(JEventLoop *loop, int eventnumber)
{
	// Get all DMCTrajectoryPoint objects for this event
	vector<const DMCTrajectoryPoint*> mctrajpoints;
	vector<const DMCThrown*> mcthrowns;
	loop->Get(mctrajpoints);
	loop->Get(mcthrowns);
	
	if(mcthrowns.size()!=1)return NOERROR;
	const DMCThrown *thrown = mcthrowns[0];
	lhit.Egen = thrown->energy();
	
	lhit.Ethrown = 0.0;
	lhit.theta_thrown = 0.0;
	lhit.E_in_bcal = 0.0;
	lhit.E_out_bcal = 0.0;
	lhit.E_out_bcal_outer = 0.0;
	lhit.E_out_bcal_inner = 0.0;
	lhit.E_leave_mother = 0.0;
	lhit.Ntraj = mctrajpoints.size();
	for(unsigned int i=0; i<mctrajpoints.size(); i++){
		const DMCTrajectoryPoint *traj = mctrajpoints[i];

		// More convenient to work with single letter variables
		double X = traj->x;
		double Y = traj->y;
		double Z = traj->z;
		double my_dE = traj->dE;
	
		// record largest energy particle's energy
		// and assume that is incident energy.
		if(traj->E>lhit.Ethrown && traj->part<=3){
			lhit.Ethrown = traj->E;
			double p = sqrt(traj->px*traj->px + traj->py*traj->py + traj->pz*traj->pz);
			lhit.theta_thrown = acos(traj->pz/p);
		}
		
		// Check if inside BCAL or not
		double R = sqrt(X*X + Y*Y);
		if(Z>=17.0 && Z<=407.0 && R>=64.3 && R<=86.172){
			// Inside BCAL
			lhit.E_in_bcal += my_dE;
		}else{
			// Outside BCAL
			lhit.E_out_bcal += my_dE;
			if(R<65.0){
				lhit.E_out_bcal_inner += my_dE;
			}else{
				lhit.E_out_bcal_outer += my_dE;
			}
		}
		
		if(fabs(X)>100.0 || fabs(Y)>=100.0 || fabs(Z)>=1000.0){
			lhit.E_leave_mother += traj->E;
		}
	}
	
	leaks->Fill();

	return NOERROR;
}

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