// $Id$
//
//    File: DL3TriggerBDT_factory.cc
// Created: Mon Aug  5 12:24:22 EDT 2013
// Creator: jrsteven (on Linux hissh0002.cmsaf.mit.edu 2.6.18-194.11.3.el5.cve20103081 x86_64)
//


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

#include <TThread.h>
#include <TStopwatch.h>
#include <TDirectory.h>
#include <TLine.h>
#include <TROOT.h>

#include "JFactoryGenerator_DL3TriggerBDT.h"
using namespace jana;

#include "DL3TriggerBDT_factory.h"

#include <TRACKING/DTrackCandidate.h>
#include <BCAL/DBCALCluster.h>
#include <BCAL/DBCALPoint.h>
#include <FCAL/DFCALCluster.h>
#include <START_COUNTER/DSCHit.h>
#include <TOF/DTOFHit.h>

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


// Initialize static members of DL3TriggerBDT_factory
static TMVA::Reader *tmvaReader = NULL;
static pthread_mutex_t tmvaReader_mutex = PTHREAD_MUTEX_INITIALIZER;

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

	l3trig = new DL3Trigger(DL3Trigger::kKEEP_EVENT, 0x0L, 0x2);
	SetFactoryFlag(PERSISTANT);
	_data.push_back(l3trig);

	// Initialize paramters
	BDT_CUT = -0.075;
	BDT_WEIGHT_PATH = "/group/halld/Users/jrsteven/2013-8-ODC/weights/";
	FILL_HISTOGRAMS = true;

	// Set parameters from command line
	gPARMS->SetDefaultParameter("L3:BDT_CUT", BDT_CUT);
	gPARMS->SetDefaultParameter("L3:BDT_WEIGHT_PATH", BDT_WEIGHT_PATH);
	gPARMS->SetDefaultParameter("L3:FILL_HISTOGRAMS", FILL_HISTOGRAMS, "Set this to 0 to disable creating and filling ROOT histograms (default is to fill them)");

	// Optionally create histograms
	if(FILL_HISTOGRAMS){
		// Initialize monitoring histograms (must be in mutex lock)
		japp->RootWriteLock();

		// Check if the first histogram has been defined already by
		// another thread. If not, then create all histograms.
		hCpuTime = (TH1F*)gROOT->FindObject("hCpuTime");
		if(!hCpuTime){
			TDirectory *main = gDirectory;
			gDirectory->mkdir("L3")->cd();

			hCpuTime = new TH1F("hCpuTime","CPU time to compute L3 decision; CPU time (ms)",100,0,100);
			hBdtClass = new TH1F("hBdtClass","Boosted Decision Tree classifier for L3 decision; BDT Classifier",100,-1.,1.);
			hTrackMomKeep = new TH1F("hTrackMomKeep","Track momentum sum for events kept by L3; Track momentum sum (GeV/c)",120,0.,12.);
			hTrackMomDiscard = new TH1F("hTrackMomDiscard","Track momentum sum for events discarded by L3; Track momentum sum (GeV/c)",120,0.,12.);

			hEbcal_EfcalKeep = new TH2F("hEbcal_EfcalKeep","BCAL energy sum vs FCAL energy sum for events kept by L3; FCAL energy sum; BCAL energy sum",120,0.,6.,120,0.,6.);
			hEbcal_EfcalDiscard = new TH2F("hEbcal_EfcalDiscard","BCAL energy sum vs FCAL energy sum for events discarded by L3; FCAL energy sum; BCAL energy sum",120,0.,6.,120,0.,6.);

			main->cd();
		}

		japp->RootUnLock();

	} // FILL_HISTOGRAMS


	// Initialize TMVA reader (do only once for all threads)
	pthread_mutex_lock(&tmvaReader_mutex);
	if(tmvaReader == NULL){
		jout<<" Initializing TMVA with the L3 trigger decision defined as a BDT response > "<<BDT_CUT<<endl;

		tmvaReader = new TMVA::Reader( "!Color:!Silent" );

		tmvaReader->AddVariable("Nstart_counter",&Nstart_counter);
		tmvaReader->AddVariable("Ntof",&Ntof);

		tmvaReader->AddVariable("Nbcal_points",&Nbcal_points);
		tmvaReader->AddVariable("Nbcal_clusters",&Nbcal_clusters);
		tmvaReader->AddVariable("EbcalPoints",&EbcalPoints);
		tmvaReader->AddVariable("EbcalClusters",&EbcalClusters);
	
		tmvaReader->AddVariable("Nfcal_clusters",&Nfcal_clusters);
		tmvaReader->AddVariable("EfcalClusters",&EfcalClusters);
	
		tmvaReader->AddVariable("Ntrack_candidates_cut",&Ntrack_candidates_cut);
		tmvaReader->AddVariable("Ptot_tracks_cut",&Ptot_tracks_cut);
   
		// book the MVA methods
		TString methodName = "BDT method";
		TString weightFile = BDT_WEIGHT_PATH; weightFile+="L3_BDT.weights.xml";
		tmvaReader->BookMVA( methodName, weightFile );
	}
	pthread_mutex_unlock(&tmvaReader_mutex);

	return NOERROR;
}

//------------------
// brun
//------------------
jerror_t DL3TriggerBDT_factory::brun(jana::JEventLoop *eventLoop, int runnumber)
{
	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);

	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t DL3TriggerBDT_factory::evnt(JEventLoop *loop, int eventnumber)
{
	TStopwatch w;
	w.Start();

	///////////////////////////////////////////////
	// Get all variables needed for L3 algorithm // 
	///////////////////////////////////////////////
	vector<const DSCHit*> schits;
	vector<const DTOFHit*> tofhits;
	vector<const DBCALPoint*> bcalpoints;
	vector<const DBCALCluster*> bcalclusters;
	vector<const DFCALCluster*> fcalclusters;
	vector<const DTrackCandidate*> candidates;

	// SC
	loop->Get(schits);
	Nstart_counter = (int)schits.size();
	
	// TOF
	loop->Get(tofhits);
	Ntof = (int)tofhits.size();

	// BCAL
	loop->Get(bcalpoints);
	Nbcal_points = (int)bcalpoints.size();
	EbcalPoints = 0.0;
	for(unsigned int i=0; i<bcalpoints.size(); i++) EbcalPoints += bcalpoints[i]->E();

	loop->Get(bcalclusters);
	Nbcal_clusters = (int)bcalclusters.size();
	EbcalClusters = 0.0;
	for(unsigned int i=0; i<bcalclusters.size(); i++) EbcalClusters += bcalclusters[i]->E();

	// FCAL
	loop->Get(fcalclusters);
	Nfcal_clusters = (int)fcalclusters.size();
	EfcalClusters = 0.0;
	for(unsigned int i=0; i<fcalclusters.size(); i++) EfcalClusters += fcalclusters[i]->getEnergy();

	// Tracks
	loop->Get(candidates);
	Ntrack_candidates_cut = 0;
	Ptot_tracks_cut = 0.0;
	TVector3 trackSum;

	for(unsigned int i=0; i<candidates.size(); i++) {
		DLorentzVector p4 = candidates[i]->lorentzMomentum();

		// Simon's POCA code for track candidates
		DVector3 doca = candidates[i]->position();
		DVector3 mom = p4.Vect();
		Float_t q = candidates[i]->charge();
		stepper->SwimToPOCAtoBeamLine(q,doca,mom);
	  
		// remove tracks with some QA criteria 
		if(candidates[i]->Ndof<3) continue;
	  
		// compute something from tracks here
		trackSum += p4.Vect();
		Ntrack_candidates_cut++;
	}
	Ptot_tracks_cut = trackSum.Mag();
	
	/////////////////////////////////
	// Assign values to DL3Trigger //
	/////////////////////////////////
	// status: should contain which criteria it failed to be rejected or that it passed (where to keep mapping?)
	// algorthim: index the current algo being used (parameter for input?)

//	DL3Trigger *l3trig = new DL3Trigger(DL3Trigger::kKEEP_EVENT, 0x0L, 0x2);
	l3trig->status      = 0x0L;
	l3trig->algorithm   = 0x2;

	// Get BDT classifier value 
	pthread_mutex_lock(&tmvaReader_mutex);
	Float_t bdtValue = tmvaReader->EvaluateMVA("BDT method");
	pthread_mutex_unlock(&tmvaReader_mutex);
	if(bdtValue > BDT_CUT)   l3trig->L3_decision = DL3Trigger::kKEEP_EVENT;
	else                     l3trig->L3_decision = DL3Trigger::kDISCARD_EVENT;

//	_data.push_back(l3trig);

	w.Stop();

	//////////////////////////////////
	// Monitoring histograms for L3 //
	//////////////////////////////////
	if(FILL_HISTOGRAMS){

		japp->RootWriteLock();
	
		hCpuTime->Fill(w.CpuTime()*1000.); //time in milli-seconds for L3 decision
		hBdtClass->Fill(bdtValue);
		if(l3trig->L3_decision == DL3Trigger::kKEEP_EVENT){
			hTrackMomKeep->Fill(Ptot_tracks_cut);
			hEbcal_EfcalKeep->Fill(EfcalClusters,EbcalClusters);
		}else {
			hTrackMomDiscard->Fill(Ptot_tracks_cut);
			hEbcal_EfcalDiscard->Fill(EfcalClusters,EbcalClusters);
		}

		japp->RootUnLock();

	} // FILL_HISTOGRAMS

	return NOERROR;
}

//------------------
// erun
//------------------
jerror_t DL3TriggerBDT_factory::erun(void)
{
	return NOERROR;
}

//------------------
// fini
//------------------
jerror_t DL3TriggerBDT_factory::fini(void)
{
	// Some additions to histograms
	TList *list;  list=hBdtClass->GetListOfFunctions();
	TLine *ln = new TLine(BDT_CUT,0,BDT_CUT,hBdtClass->GetMaximum()); ln->SetLineColor(kRed); list->Add(ln);

	// Delete TMVA reader
	pthread_mutex_lock(&tmvaReader_mutex);
	if(tmvaReader){
		delete tmvaReader;
		tmvaReader = NULL;
	}
	pthread_mutex_unlock(&tmvaReader_mutex);
	
	return NOERROR;
}