#include <TRACKING/DMCThrown.h>

#include <FDC/DFDCHit.h>
#include <FDC/DFDCPseudo.h>
#include <CDC/DCDCHit.h>
#include <TRACKING/DTrackCandidate.h>
#include <TRACKING/DTrackWireBased.h>
#include <TRACKING/DTrackTimeBased.h>
#include <PID/DChargedTrack.h>
#include <PID/DChargedTrackHypothesis.h>
#include <BCAL/DBCALHit.h>
#include <BCAL/DBCALUnifiedHit.h>
#include <BCAL/DBCALCluster.h>
#include <FCAL/DFCALHit.h>
#include <FCAL/DFCALCluster.h>

#include <START_COUNTER/DSCHit.h>
#include <TOF/DTOFHit.h>
#include <TOF/DTOFPoint.h>
#include <BCAL/DBCALShower.h>
#include <FCAL/DFCALShower.h>

#include <PID/DBeamPhoton.h>
#include <TRIGGER/DMCTrigger.h>

#include <TAGGER/DTAGMHit.h>
#include <TAGGER/DTAGHHit.h>
#include <PID/DEventRFBunch.h>

#include "DEventProcessor_L3.h"
#include "DFactoryGenerator_EventRFBunch_L3.h"

// The executable should define the ROOTfile global variable. It will
// be automatically linked when dlopen is called.
extern TFile *ROOTfile;

// Routine used to create our DEventProcessor
extern "C"
{
        void InitPlugin(JApplication *app)
        {
                InitJANAPlugin(app);
                app->AddProcessor(new DEventProcessor_L3());
		app->AddFactoryGenerator(new DFactoryGenerator_EventRFBunch_L3());
        }
} // "C"

//define static local variable //declared in header file
thread_local DTreeFillData DEventProcessor_L3::dTreeFillData;

//------------------
// DEventProcessor_L3 (Constructor)
//------------------
DEventProcessor_L3::DEventProcessor_L3()
{

}

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

}

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

	// Create histograms here
	hMicroDeltaTAll = new TH2F("hMicroDeltaTAll","hMicroDeltaTAll; Microscope E_{#gamma} (GeV); #Deltat_{Micro #gamma - thrown RF} (ns)",120,0.,12.,100,-10,10);
	hHodoDeltaTAll = new TH2F("hHodoDeltaTAll","hHodoDeltaTAll; Hodoscope E_{#gamma} (GeV); #Deltat_{Hodo #gamma - thrown RF} (ns)",120,0.,12.,100,-10,10);

	hMicroTrueRFTag = new TH1F("hMicroTrueRFTag","hMicroTrueRFTag; # Microscope Tags",20,0,20);
	hMicroWrongRFTag = new TH1F("hMicroWrongRFTag","hMicroWrongRFTag; # Microscope Tags",20,0,20);
	hHodoTrueRFTag = new TH1F("hHodoTrueRFTag","hHodoTrueRFTag; # Hodoscope Tags",20,0,20);
	hHodoWrongRFTag = new TH1F("hHodoWrongRFTag","hHodoWrongRFTag; # Hodoscope Tags",20,0,20);
	hHodoTrueRFTagHigh = new TH1F("hHodoTrueRFTagHigh","hHodoTrueRFTagHigh; # Hodoscope Tags High",20,0,20);
	hHodoWrongRFTagHigh = new TH1F("hHodoWrongRFTagHigh","hHodoWrongRFTagHigh; # Hodoscope Tags High",20,0,20);
	hCombinedTrueRFTag = new TH1F("hCombinedTrueRFTag","hCombinedTrueRFTag; # Combined Tags",20,0,20);
	hCombinedWrongRFTag = new TH1F("hCombinedWrongRFTag","hCombinedWrongRFTag; # Combined Tags",20,0,20);

	hDeltaT = new TH1F("hDeltaT","hDeltaT; #Deltat_{Beam #gamma - RF} (ns)",100,-10,10);
	hDeltaT_Egamma = new TH2F("hDeltaT_Egamma","hDeltaT_Egamma; E_{#gamma} (GeV); #Deltat_{Beam #gamma - RF} (ns)",120,0.,12.,100,-10,10);
	hDeltaT_Type = new TH2F("hDeltaT_Type","hDeltaT_Type; RFBunch ID Type; #Deltat_{Beam #gamma - RF} (ns)",5,0.,5.,100,-10,10);
	hEgamma_Type = new TH2F("hEgamma_Type","hEgamma_Type;  RFBunch ID Type; E_{#gamma} (GeV)",5,0.,5.,120,0.,12.);

	//TTREE INTERFACE
        //MUST DELETE WHEN FINISHED: OR ELSE DATA WON'T BE SAVED!!!
        dTreeInterface = DTreeInterface::Create_DTreeInterface("TrainingL3_Tree", "tree_TrainingL3.root");

        //TTREE BRANCHES
        DTreeBranchRegister locTreeBranchRegister;

	// MC variables
        locTreeBranchRegister.Register_Single<float>("Ephoton_truth");
	Nparticle_types = 50;
	locTreeBranchRegister.Register_Single<int>("Nparticle_types");
	locTreeBranchRegister.Register_FundamentalArray<int>("numParticlesGeneratedByType", "Nparticle_types", Nparticle_types);

	int locNumTrackThrown = 100;
	locTreeBranchRegister.Register_Single<int>("Ntrack_thrown");
	locTreeBranchRegister.Register_FundamentalArray<int>("Track_thrown_PID", "Ntrack_thrown", locNumTrackThrown);
	locTreeBranchRegister.Register_FundamentalArray<int>("Track_thrown_charge", "Ntrack_thrown", locNumTrackThrown);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_thrown_px", "Ntrack_thrown", locNumTrackThrown);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_thrown_py", "Ntrack_thrown", locNumTrackThrown);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_thrown_pz", "Ntrack_thrown", locNumTrackThrown);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_thrown_E", "Ntrack_thrown", locNumTrackThrown);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_thrown_doca_x", "Ntrack_thrown", locNumTrackThrown);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_thrown_doca_y", "Ntrack_thrown", locNumTrackThrown);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_thrown_doca_z", "Ntrack_thrown", locNumTrackThrown);

	// Reco variables
	locTreeBranchRegister.Register_Single<bool>("L1_fired");
	locTreeBranchRegister.Register_Single<float>("Nstart_counter");
	locTreeBranchRegister.Register_Single<float>("Ntof");
	locTreeBranchRegister.Register_Single<float>("Ntof_point");

	int locNumBCALshowers = 100;
	locTreeBranchRegister.Register_Single<int>("Nbcal_showers");
	locTreeBranchRegister.Register_Single<float>("EbcalShowers");
	locTreeBranchRegister.Register_FundamentalArray<double>("BCAL_shower_x", "Nbcal_showers", locNumBCALshowers);
	locTreeBranchRegister.Register_FundamentalArray<double>("BCAL_shower_y", "Nbcal_showers", locNumBCALshowers);
	locTreeBranchRegister.Register_FundamentalArray<double>("BCAL_shower_z", "Nbcal_showers", locNumBCALshowers);
	locTreeBranchRegister.Register_FundamentalArray<double>("BCAL_shower_E", "Nbcal_showers", locNumBCALshowers);
	locTreeBranchRegister.Register_FundamentalArray<double>("BCAL_shower_t0", "Nbcal_showers", locNumBCALshowers);

	int locNumFCALshowers = 100;
	locTreeBranchRegister.Register_Single<int>("Nfcal_showers");
	locTreeBranchRegister.Register_Single<float>("EfcalShowers");
	locTreeBranchRegister.Register_FundamentalArray<double>("FCAL_shower_x", "Nfcal_showers", locNumFCALshowers);
	locTreeBranchRegister.Register_FundamentalArray<double>("FCAL_shower_y", "Nfcal_showers", locNumFCALshowers);
	locTreeBranchRegister.Register_FundamentalArray<double>("FCAL_shower_z", "Nfcal_showers", locNumFCALshowers);
	locTreeBranchRegister.Register_FundamentalArray<double>("FCAL_shower_E", "Nfcal_showers", locNumFCALshowers);
	locTreeBranchRegister.Register_FundamentalArray<double>("FCAL_shower_t0", "Nfcal_showers", locNumFCALshowers);

	locTreeBranchRegister.Register_Single<float>("NtagRF");
	locTreeBranchRegister.Register_Single<float>("Ntagm");
	locTreeBranchRegister.Register_Single<float>("Ntagh");

	int locNumTrackCandidates = 100;
	locTreeBranchRegister.Register_Single<int>("Ntrack_candidates");
	locTreeBranchRegister.Register_Single<float>("Ptot_tracks_cand");
	locTreeBranchRegister.Register_FundamentalArray<int>("Track_cand_charge", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<int>("Track_cand_NDF", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_cand_ChiSq", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_cand_px", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_cand_py", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_cand_pz", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_cand_E", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_cand_doca_x", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_cand_doca_y", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_cand_doca_z", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<int>("Track_candWire_charge", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<int>("Track_candWire_NDF", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_candWire_ChiSq", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_candWire_px", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_candWire_py", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_candWire_pz", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_candWire_E", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_candThrown_px", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_candThrown_py", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_candThrown_pz", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_candThrown_E", "Ntrack_candidates", locNumTrackCandidates);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_candThrown_mass", "Ntrack_candidates", locNumTrackCandidates);

	int locNumTrackWires = 100;
	locTreeBranchRegister.Register_Single<int>("Ntrack_wires");
	locTreeBranchRegister.Register_Single<float>("Ptot_tracks_wire");
	locTreeBranchRegister.Register_FundamentalArray<int>("Track_wire_charge", "Ntrack_wires", locNumTrackWires);
	locTreeBranchRegister.Register_FundamentalArray<int>("Track_wire_NDF", "Ntrack_wires", locNumTrackWires);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_wire_ChiSq", "Ntrack_wires", locNumTrackWires);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_wire_px", "Ntrack_wires", locNumTrackWires);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_wire_py", "Ntrack_wires", locNumTrackWires);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_wire_pz", "Ntrack_wires", locNumTrackWires);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_wire_E", "Ntrack_wires", locNumTrackWires);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_wire_doca_x", "Ntrack_wires", locNumTrackWires);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_wire_doca_y", "Ntrack_wires", locNumTrackWires);
	locTreeBranchRegister.Register_FundamentalArray<double>("Track_wire_doca_z", "Ntrack_wires", locNumTrackWires);
	
	//REGISTER BRANCHES
        dTreeInterface->Create_Branches(locTreeBranchRegister);
	
	// tracking code here
	swim_steps = new DReferenceTrajectory::swim_step_t[max_swim_steps];

	return NOERROR;
}

//------------------
// brun
//------------------
jerror_t DEventProcessor_L3::brun(JEventLoop *eventLoop, int32_t runnumber)
{
        // Get attenuation parameters
        double L_over_2 = DBCALGeometry::GetBCAL_length()/2.0;
        double Xo = DBCALGeometry::ATTEN_LENGTH;
        unattenuate_to_center = exp(+L_over_2/Xo);

	dApplication = dynamic_cast<DApplication*>(eventLoop->GetJApplication());

	// Special magnetic field stepper for POCA to beamline
	bfield = dApplication->GetBfield();
	stepper = new DMagneticFieldStepper(bfield);
	stepper->SetStepSize(1.0);
       
	// Get fitters
	vector <const DTrackFitter*> fitters;
	eventLoop->Get(fitters);
	fitter = const_cast<DTrackFitter*>(fitters[0]);

	// reference trajectory for fitter
	rt = new DReferenceTrajectory(bfield, 1.0, swim_steps, max_swim_steps);

	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t DEventProcessor_L3::evnt(JEventLoop *loop, uint64_t eventnumber)
{
	// get MC truth info
	getMCTruth(loop);

	// start with fast detectors for L3
	vector<const DSCHit*> schits;
	vector<const DTOFHit*> tofhits;
	vector<const DTOFPoint*> tofpoints;
	vector<const DBCALShower*> bcalshowers;
	vector<const DFCALShower*> fcalshowers;

	// SC
	loop->Get(schits);
	int Nstart_counter = (int)schits.size();
	dTreeFillData.Fill_Single<float>("Nstart_counter", Nstart_counter);

	// TOF
	loop->Get(tofhits);
	int Ntof = (int)tofhits.size();
	dTreeFillData.Fill_Single<float>("Ntof", Ntof);
	loop->Get(tofpoints);
	int Ntof_point = (int)tofpoints.size();
	dTreeFillData.Fill_Single<float>("Ntof_point", Ntof_point);

	// BCAL showers
	loop->Get(bcalshowers);
	int Nbcal_showers = (int)bcalshowers.size();
	dTreeFillData.Fill_Single<int>("Nbcal_showers", Nbcal_showers);
	double EbcalShowers = 0.0;
	for(unsigned int i=0; i<bcalshowers.size(); i++) {
		EbcalShowers += bcalshowers[i]->E;
		dTreeFillData.Fill_Array<double>("BCAL_shower_x", bcalshowers[i]->x, i);
		dTreeFillData.Fill_Array<double>("BCAL_shower_y", bcalshowers[i]->y, i);
		dTreeFillData.Fill_Array<double>("BCAL_shower_z", bcalshowers[i]->z, i);
		dTreeFillData.Fill_Array<double>("BCAL_shower_E", bcalshowers[i]->E, i);
		dTreeFillData.Fill_Array<double>("BCAL_shower_t0", bcalshowers[i]->t, i);
	}
	dTreeFillData.Fill_Single<float>("EbcalShowers", EbcalShowers);

	// FCAL showers
        loop->Get(fcalshowers);
        int Nfcal_showers = (int)fcalshowers.size();
	dTreeFillData.Fill_Single<int>("Nfcal_showers", Nfcal_showers);
        double EfcalShowers = 0.0;
        for(unsigned int i=0; i<fcalshowers.size(); i++) {
		EfcalShowers += fcalshowers[i]->getEnergy();
		dTreeFillData.Fill_Array<double>("FCAL_shower_x", fcalshowers[i]->getPosition().X(), i);
		dTreeFillData.Fill_Array<double>("FCAL_shower_y", fcalshowers[i]->getPosition().Y(), i);
		dTreeFillData.Fill_Array<double>("FCAL_shower_z", fcalshowers[i]->getPosition().Z(), i);
		dTreeFillData.Fill_Array<double>("FCAL_shower_E", fcalshowers[i]->getEnergy(), i);
		dTreeFillData.Fill_Array<double>("FCAL_shower_t0", fcalshowers[i]->getTime(), i);
	}
	dTreeFillData.Fill_Single<float>("EfcalShowers", EfcalShowers);

	// get L1 trigger
	bool L1_fired = getL1(loop);

	// get tracks info
	getTracks(loop);

	// get tagger info
	getTagger(loop, L1_fired);

	// write L3 trigger tree
	dTreeInterface->Fill(dTreeFillData);

	return NOERROR;
}

//------------------
bool DEventProcessor_L3::getL1(JEventLoop *loop)
{

	////////////////////////
	// copy of L1 trigger //
	////////////////////////

	vector<const DBCALUnifiedHit*> bcalhits_unified;
	vector<const DFCALHit*> fcalhits;
	loop->Get(bcalhits_unified);
	loop->Get(fcalhits);

	double BCAL_Eupstream = 0.0;
        double BCAL_Edownstream = 0.0;
	for(unsigned int i=0; i< bcalhits_unified.size(); i++){
                const DBCALUnifiedHit* bcalhit_unified = bcalhits_unified[i];
		if(bcalhit_unified->E < 0.020)   //removes dark pulse effects (new in mcsmear) but lowers efficiency
			continue;

                if(bcalhit_unified->end == DBCALGeometry::kUpstream){
                        BCAL_Eupstream += bcalhit_unified->E;
                }else{
                        BCAL_Edownstream += bcalhit_unified->E;
                }
        }
	
        // Calculate "energy" sum for BCAL in GeV-ish units
        double EbcalTrig = unattenuate_to_center * (BCAL_Eupstream + BCAL_Edownstream)/2.0;

	// FCAL
        loop->Get(fcalhits);
        double Efcal = 0.0;
        for(unsigned int i=0; i<fcalhits.size(); i++) Efcal += fcalhits[i]->E;

        // BCAL and FCAL
	bool L1_fired = false;
        bool sum_cut = (EbcalTrig + 4.0*Efcal)>=2.0;
        if(sum_cut && EbcalTrig>0.200 && Efcal>0.030) 
		L1_fired = true;

	dTreeFillData.Fill_Single<bool>("L1_fired", L1_fired);

	return L1_fired;
}

//------------------
jerror_t DEventProcessor_L3::getTracks(JEventLoop *loop)
{

	///////////////////////////
	// tracking needs for L3 //
	///////////////////////////
	vector<const DTrackCandidate*> candidates;

	// Tracks
	loop->Get(candidates);
	int Ntrack_candidates = (int)candidates.size();
	dTreeFillData.Fill_Single<int>("Ntrack_candidates", Ntrack_candidates);

	double Ptot_tracks_cand = 0.0;
	for(unsigned int i=0; i<candidates.size(); i++) {
		Ptot_tracks_cand += candidates[i]->pmag();
		dTreeFillData.Fill_Array<int>("Track_cand_charge", (int)candidates[i]->charge(), i);
		dTreeFillData.Fill_Array<int>("Track_cand_NDF", candidates[i]->Ndof, i);
		dTreeFillData.Fill_Array<double>("Track_cand_ChiSq", candidates[i]->chisq, i);
		dTreeFillData.Fill_Array<double>("Track_cand_px", candidates[i]->lorentzMomentum().X(), i);
		dTreeFillData.Fill_Array<double>("Track_cand_py", candidates[i]->lorentzMomentum().Y(), i);
		dTreeFillData.Fill_Array<double>("Track_cand_pz", candidates[i]->lorentzMomentum().Z(), i);
		dTreeFillData.Fill_Array<double>("Track_cand_E", candidates[i]->lorentzMomentum().E(), i);

/*
		// do wire-based track refit
		float mass = 0.1396;
		fitter->SetFitType(DTrackFitter::kWireBased);
		// Swim a reference trajectory using the candidate starting momentum and position
		rt->SetMass(mass);
		rt->Swim(candidates[i]->position(),candidates[i]->momentum(),candidates[i]->charge());
		DTrackFitter::fit_status_t status = fitter->FindHitsAndFitTrack(*candidates[i],rt,loop,mass,candidates[i]->Ndof+3);
		
		// wire-based refit failed
		if(status == DTrackFitter::kFitNotDone || status == DTrackFitter::kFitFailed){
			dTreeFillData.Fill_Array<int>("Track_candWire_NDF", -999, i);
			//cout<<"bad fit routine"<<endl;
		}
		else{ // wire-based refit succeeeded
			DKinematicData kd = fitter->GetFitParameters();
			dTreeFillData.Fill_Array<int>("Track_candWire_charge", (int)kd.charge(), i);
			dTreeFillData.Fill_Array<int>("Track_candWire_NDF", fitter->GetNdof(), i);
			dTreeFillData.Fill_Array<double>("Track_candWire_ChiSq", fitter->GetChisq(), i);
			dTreeFillData.Fill_Array<double>("Track_candWire_px", kd.lorentzMomentum().X(), i);
			dTreeFillData.Fill_Array<double>("Track_candWire_py", kd.lorentzMomentum().X(), i);
			dTreeFillData.Fill_Array<double>("Track_candWire_pz", kd.lorentzMomentum().X(), i);
			dTreeFillData.Fill_Array<double>("Track_candWire_E", kd.lorentzMomentum().E(), i);
			//cout<<" candidate momentum  = "<<candidates[i]->momentum().Mag()<<endl;
			//cout<<" wire-based momentum = "<<kd.momentum().Mag()<<endl;
		}
*/

		// Simon's POCA code
		DVector3 pos = candidates[i]->position();
		DVector3 mom = candidates[i]->lorentzMomentum().Vect();
		Float_t q = candidates[i]->charge();
		stepper->SwimToPOCAtoBeamLine(q,pos,mom);
		dTreeFillData.Fill_Array<double>("Track_cand_doca_x", pos.X(), i);
		dTreeFillData.Fill_Array<double>("Track_cand_doca_y", pos.Y(), i);
		dTreeFillData.Fill_Array<double>("Track_cand_doca_z", pos.Z(), i);
			
		vector<const DMCThrown*> locMCThrownsFinalState;
		loop->Get(locMCThrownsFinalState, "FinalState");
		
		// Match to thrown tracks
		float maxFOM = -1;
		for(size_t loc_i = 0; loc_i < locMCThrownsFinalState.size(); ++loc_i){
			const DMCThrown *locMCThrown = locMCThrownsFinalState[loc_i];
			
			//only get write sign tracks
			if(candidates[i]->charge() != ParticleCharge((Particle_t)(locMCThrown->type))) continue; 
			
			float MatchFOM = Calc_MatchFOM(locMCThrown->momentum(), candidates[i]->momentum());
			
			if(MatchFOM > maxFOM){
				dTreeFillData.Fill_Array<double>("Track_candThrown_px", locMCThrown->momentum().X(), i);
				dTreeFillData.Fill_Array<double>("Track_candThrown_py", locMCThrown->momentum().Y(), i);
				dTreeFillData.Fill_Array<double>("Track_candThrown_pz", locMCThrown->momentum().Z(), i);
				dTreeFillData.Fill_Array<double>("Track_candThrown_E", locMCThrown->energy(), i);
				dTreeFillData.Fill_Array<double>("Track_candThrown_mass", locMCThrown->mass(), i);
			}
		}
	}       
	dTreeFillData.Fill_Single<float>("Ptot_tracks_cand", Ptot_tracks_cand);
	
	// wire based tracking
	vector<const DTrackTimeBased*> tracks_wb;
	loop->Get(tracks_wb);
	int Ntrack_wires = (int)tracks_wb.size();
	dTreeFillData.Fill_Single<int>("Ntrack_wires", Ntrack_wires);

	double Ptot_tracks_wire = 0.;
	for(unsigned int i=0; i<tracks_wb.size(); i++) {
		Ptot_tracks_wire += tracks_wb[i]->pmag();
		dTreeFillData.Fill_Array<int>("Track_wire_charge", (int)tracks_wb[i]->charge(), i);
		dTreeFillData.Fill_Array<int>("Track_wire_NDF", tracks_wb[i]->Ndof, i);
		dTreeFillData.Fill_Array<double>("Track_wire_ChiSq", tracks_wb[i]->chisq, i);
		dTreeFillData.Fill_Array<double>("Track_wire_px", tracks_wb[i]->lorentzMomentum().X(), i);
		dTreeFillData.Fill_Array<double>("Track_wire_py", tracks_wb[i]->lorentzMomentum().Y(), i);
		dTreeFillData.Fill_Array<double>("Track_wire_pz", tracks_wb[i]->lorentzMomentum().Z(), i);
		dTreeFillData.Fill_Array<double>("Track_wire_E", tracks_wb[i]->lorentzMomentum().E(), i);

		// Simon's POCA code
		DVector3 pos = tracks_wb[i]->position();
		DVector3 mom = tracks_wb[i]->lorentzMomentum().Vect();
		Float_t q = tracks_wb[i]->charge();
		stepper->SwimToPOCAtoBeamLine(q,pos,mom);
		dTreeFillData.Fill_Array<double>("Track_wire_doca_x", pos.X(), i);
		dTreeFillData.Fill_Array<double>("Track_wire_doca_y", pos.Y(), i);
		dTreeFillData.Fill_Array<double>("Track_wire_doca_z", pos.Z(), i);
	}
	dTreeFillData.Fill_Single<float>("Ptot_tracks_wire", Ptot_tracks_wire);

	return NOERROR;
}


//------------------
jerror_t DEventProcessor_L3::getTagger(JEventLoop *loop, bool L1_fired)
{

	/////////////////////////
	// tagger needs for L3 //
	/////////////////////////

	vector<const DBeamPhoton*> beamphotons;
	loop->Get(beamphotons, "MCGEN");
	double Ephoton_truth = beamphotons.size()>0 ? beamphotons[0]->energy():0.0;

	// microscope and hodoscope hits
	vector<const DTAGMHit*> microHits;
	vector<const DTAGHHit*> hodoHits;
	loop->Get(microHits);
	loop->Get(hodoHits);
	int nMicroTrueRF = 0;
	int nMicroWrongRF = 0;
	int nHodoTrueRF = 0;
	int nHodoWrongRF = 0;
	int nHodoTrueRFHigh = 0;
	int nHodoWrongRFHigh = 0;
	
	// true RF bunch for comparison
	vector<const DEventRFBunch*> eventRFBunchThrown;
	loop->Get(eventRFBunchThrown, "Thrown");

	// time window to select RF bunch
	float lowTrue = -2.0;
	float highTrue = 2.0;
	float lowWrong = 2.0;
	float highWrong = 6.0;

	// energy window to select micro and hodoscope (high energy)
	float lowMicro = 8.2;
	float highMicro = 9.2;
	float lowHodo = 9.2;
	float highHodo = 12;

	// FILL HISTOGRAMS
	// Since we are filling histograms local to this plugin, it will not interfere with other ROOT operations: can use plugin-wide ROOT fill lock
	japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK

	// fill tagger histograms for L1 triggered events (counts true and accidentals)
	if(L1_fired) {
		for(size_t loc_i = 0; loc_i < microHits.size(); ++loc_i){
			double locDeltaT = microHits[loc_i]->t - eventRFBunchThrown[0]->dTime;
			hMicroDeltaTAll->Fill(microHits[loc_i]->E,locDeltaT);
			if(locDeltaT > lowTrue && locDeltaT < highTrue) // pick true RF bucket from MC event
				nMicroTrueRF++; 
			if(locDeltaT > lowWrong && locDeltaT < highWrong) // pick wrong RF bucket not from MC event 
				nMicroWrongRF++;
		}
		if(Ephoton_truth > lowMicro && Ephoton_truth < highMicro) 
			hMicroTrueRFTag->Fill(nMicroTrueRF);
		hMicroWrongRFTag->Fill(nMicroWrongRF);
		
		for(size_t loc_i = 0; loc_i < hodoHits.size(); ++loc_i){
			double locDeltaT = hodoHits[loc_i]->t - eventRFBunchThrown[0]->dTime;
			hHodoDeltaTAll->Fill(hodoHits[loc_i]->E,locDeltaT);
			if(locDeltaT > lowTrue && locDeltaT < highTrue) { // pick true RF bucket from MC event
				nHodoTrueRF++; 
				if(hodoHits[loc_i]->E > lowHodo) 
					nHodoTrueRFHigh++;
			}	
			if(locDeltaT > lowWrong && locDeltaT < highWrong) { // pick RF bucket not from true MC event
				nHodoWrongRF++; 
				if(hodoHits[loc_i]->E > lowHodo) 
					nHodoWrongRFHigh++;
			}
		}
		if(Ephoton_truth > 3.0 && Ephoton_truth < highHodo && !(Ephoton_truth > lowMicro && Ephoton_truth < highMicro))
			hHodoTrueRFTag->Fill(nHodoTrueRF);
		if(Ephoton_truth > lowHodo && Ephoton_truth < highHodo)
			hHodoTrueRFTagHigh->Fill(nHodoTrueRFHigh);
		hHodoWrongRFTag->Fill(nHodoWrongRF);
		hHodoWrongRFTagHigh->Fill(nHodoWrongRFHigh);
		
		if(Ephoton_truth > lowMicro && Ephoton_truth < highHodo)
			hCombinedTrueRFTag->Fill(nMicroTrueRF+nHodoTrueRFHigh);
		hCombinedWrongRFTag->Fill(nMicroWrongRF+nHodoWrongRFHigh);
	}
	japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK

	// candidate beam photons (same as sum of microscope and hodoscope hits)
	vector<const DBeamPhoton*> beamphotonCandidates;
	loop->Get(beamphotonCandidates);

	// select candidate RF bunches (not using any tracking)
	vector<const DEventRFBunch*> eventRFBunch;
	loop->Get(eventRFBunch, "L3");
	int type = 0;
	if(eventRFBunch[0]->dNumParticleVotes >= 1) type = 1;

	// monitor how well RF bunch is selected for L1 triggered events
	japp->RootFillLock(this); //ACQUIRE ROOT FILL LOCK
	double locDeltaT = eventRFBunchThrown[0]->dTime - eventRFBunch[0]->dTime;
	if(L1_fired) {
		hEgamma_Type->Fill(type,Ephoton_truth);
		if(Ephoton_truth > lowMicro) hDeltaT_Type->Fill(type,locDeltaT);
		if(eventRFBunch[0]->dNumParticleVotes >= 1){
			hDeltaT->Fill(locDeltaT);
			hDeltaT_Egamma->Fill(Ephoton_truth,locDeltaT);
		}
	}
	japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK

	// store variables to use in evaluating L3 trigger (not used in BDT training)
	int tagRF = 0;
	int hitMicroscope=0;
	int hitHodoscope=0;
	if(eventRFBunch[0]->dNumParticleVotes >= 1){ // only cut on tagger if can find RF bunch with high energy calorimeter cluster
		tagRF = 1;
		
		for(size_t loc_k = 0; loc_k < beamphotonCandidates.size(); ++loc_k){
			double locDeltaT = beamphotonCandidates[loc_k]->time() - eventRFBunch[0]->dTime;
			
			// only look in 2 RF bunches which contain 99% of true beam photons 
			if(locDeltaT > 2.004 || locDeltaT < -2.004) 
				continue;
			
			const DBeamPhoton* beamphotonCandidate = beamphotonCandidates[loc_k];

			if(beamphotonCandidate->energy() > lowMicro && beamphotonCandidate->energy() < highMicro) 
				hitMicroscope++;
			else if(beamphotonCandidate->energy() > lowHodo && beamphotonCandidate->energy() < highHodo) 
				hitHodoscope++;
		}
	}
	dTreeFillData.Fill_Single<float>("NtagRF", tagRF);
	dTreeFillData.Fill_Single<float>("Ntagm", hitMicroscope);
	dTreeFillData.Fill_Single<float>("Ntagh", hitHodoscope);

	return NOERROR;
}

//------------------
// erun
//------------------
jerror_t DEventProcessor_L3::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_L3::fini(void)
{
	// If another DEventProcessor is in the list ahead of this one, then
	// it will have finished before this is called. e.g. closed the
	// ROOT file!

        delete[] swim_steps;

	delete dTreeInterface; //saves trees to file, closes file

	return NOERROR;
}

//------------------
jerror_t DEventProcessor_L3::getMCTruth(JEventLoop *loop)
{
	// beam photons
	vector<const DBeamPhoton*> beamphotons;
	loop->Get(beamphotons, "MCGEN");
	double Ephoton_truth = beamphotons.size()>0 ? beamphotons[0]->energy():0.0;
	dTreeFillData.Fill_Single<float>("Ephoton_truth", Ephoton_truth);

	vector<const DMCThrown*> locMCThrowns;
        loop->Get(locMCThrowns);

	vector<const DMCThrown*> locMCThrownsFinalState;
        loop->Get(locMCThrownsFinalState, "FinalState");

	Int_t numParticlesGeneratedByType[Nparticle_types];
	for(unsigned int i=0;i<Nparticle_types;i++){
		numParticlesGeneratedByType[i] = 0;
	}
	
	int Ntrack_thrown=0;
	for(size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i)
        {
		const DMCThrown *thrown = locMCThrowns[loc_i];
		
		//skip muons, electrons and neutrinos (hack for now)
		if(thrown->type >= 2 && thrown->type <=6) continue;
		
		//skip bad thrown particles from geant interactions?
		if(thrown->type == 14 || thrown->type == 15 || thrown->type == 8 || thrown->type == 9 || thrown->type == 11 || thrown->type == 12){
			if(thrown->myid>=2000 && (thrown->position().Perp()>10. || thrown->position().Z()<0. || thrown->position().Z()>200)) continue; 
			if(thrown->myid>=3000) continue;
		}
		
		// Save how many particles of each ParticleType were created
		int type = thrown->type;
		numParticlesGeneratedByType[(Int_t)type]++;

		dTreeFillData.Fill_Array<int>("Track_thrown_PID", (int)thrown->PID(), Ntrack_thrown);
		dTreeFillData.Fill_Array<int>("Track_thrown_charge", (int)thrown->charge(), Ntrack_thrown);
		dTreeFillData.Fill_Array<double>("Track_thrown_px", thrown->lorentzMomentum().X(), Ntrack_thrown);
		dTreeFillData.Fill_Array<double>("Track_thrown_py", thrown->lorentzMomentum().Y(), Ntrack_thrown);
		dTreeFillData.Fill_Array<double>("Track_thrown_pz", thrown->lorentzMomentum().Z(), Ntrack_thrown);
		dTreeFillData.Fill_Array<double>("Track_thrown_E", thrown->lorentzMomentum().E(), Ntrack_thrown);
		dTreeFillData.Fill_Array<double>("Track_thrown_doca_x", thrown->position().X(), Ntrack_thrown);
		dTreeFillData.Fill_Array<double>("Track_thrown_doca_y", thrown->position().Y(), Ntrack_thrown);
		dTreeFillData.Fill_Array<double>("Track_thrown_doca_z", thrown->position().Z(), Ntrack_thrown);
		Ntrack_thrown++;
	}
	dTreeFillData.Fill_Single<int>("Nparticle_types", Nparticle_types);
	dTreeFillData.Fill_Single<int>("Ntrack_thrown", Ntrack_thrown);
	
	for(unsigned int i=0;i<Nparticle_types;i++){
		dTreeFillData.Fill_Array<int>("numParticlesGeneratedByType", numParticlesGeneratedByType[i], i);
	}
	
	return NOERROR;
}

double DEventProcessor_L3::Calc_MatchFOM(const DVector3& locMomentum_Thrown, const DVector3& locMomentum_Detected) 
{
        double locDeltaPOverP = (locMomentum_Thrown.Mag() > 0.0) ? (locMomentum_Thrown.Mag() - locMomentum_Detected.Mag())/locMomentum_Thrown.Mag() : 9.9E9; //mean is zero
        double locDeltaPOverPSigma = 0.03;
        double locDeltaPOverPChiSq = locDeltaPOverP*locDeltaPOverP/(locDeltaPOverPSigma*locDeltaPOverPSigma);
	
        double locDeltaTheta = locMomentum_Thrown.Angle(locMomentum_Detected); //mean is zero
        double locDeltaThetaSigma = 0.5*TMath::Pi()/180.0; //0.5 degrees: phi dominates (theta is ~0.03 degrees)
        double locDeltaThetaChiSq = locDeltaTheta*locDeltaTheta/(locDeltaThetaSigma*locDeltaThetaSigma);

        double locFOM = TMath::Prob(locDeltaPOverPChiSq + locDeltaThetaChiSq, 2);

        return locFOM;
}