// $Id$
//
//    File: MyProcessor.cc
// Created: Mon Aug 28 EDT 2006
// Creator: davidl (on Darwin swire-b241.jlab.org 8.7.0 powerpc)
//

#include <sstream>
using namespace std;

#include <TThread.h>

#include <JANA/JFactory.h>

#include "DEventSourceET/DADC.h"
#include "DEventSourceET/DTDC.h"
#include "DEventSourceET/DTrigger.h"
#include "DEventSourceET/DEPICS.h"
#include "BCAL/DBCAL_ADCHit_factory.h"
#include "BCAL/DBCAL_TDCHit.h"
#include "TAGGER/DTAGGER_ETCoinc_factory.h"
#include "TRIGGER/DTRIGGER_TDCHit.h"

#include "MyProcessor.h"

// This program is used to create ROOT Tree DSTs of the 
// bcal06 data.

// This global is set to the mean time (timewalk corrected)
// of BCAL cell 8. It is used in TAGGERTDCsort below to
// order tagger "photons" such that the most in-time with
// the BCAL is first.
double BCAL_TIME; 

// Sorting functions used to put TDC hits in order
bool TrigTDCsort(const DTRIGGER_TDCHit* const &t1, const DTRIGGER_TDCHit* const &t2) {
	if(t1->t != t2->t)return fabs(t1->t) < fabs(t2->t);
	return t1->tid < t1->tid;
}

// Sorting functions used to put TDC hits in order
bool TAGGERTDCsort(const DTAGGER_ETCoinc* const &t1, const DTAGGER_ETCoinc* const &t2) {
	if(t1->time != t2->time){
		return fabs(t1->time+BCAL_TIME-531.0) < fabs(t2->time+BCAL_TIME-531.0);
	}
	return t1->T_id < t1->T_id;
}


//------------------
// init
//------------------
MyProcessor::MyProcessor(void)
{
	// ROOT file is opened in brun so we can use runnumber 
	// in the filename. If you want one file for several
	// runs (e.g. multiple files from different runs are
	// specified on the command line) then open the file
	// here.
	file = NULL;
	runnumber = 0;
	pthread_mutex_init(&root_mutex,NULL);
}

//------------------
// brun
//------------------
jerror_t MyProcessor::brun(JEventLoop *eventLoop, int runnumber)
{
	// Special events (like prestart and go) will not have run number 
	// info in them. Ignore these.
	if(runnumber==0 || runnumber==this->runnumber)return NOERROR;

	// Lock access to ROOT
	pthread_mutex_lock(&root_mutex);

	// Create ROOT file
	char fname[256];
	sprintf(fname,"bcal_dst%05d.root", runnumber);
	file = new TFile(fname, "RECREATE");
	cout<<"Opened \""<<fname<<"\" ...."<<endl;
	this->runnumber = runnumber;

	// Create TTrees
	bcal = new TTree("bcal", "BCAL raw values",8*1024*1024);
	stringstream sN, sS, stN, stS, stmean;
	for(int i=1; i<=18; i++){
		sN<<"n"<<i<<"/S";
		sS<<"s"<<i<<"/S";
		stN<<"tn"<<i<<"/F";
		stS<<"ts"<<i<<"/F";
		stmean<<"tmean"<<i<<"/F";
		if(i!=18){
			sN<<":";
			sS<<":";
			stN<<":";
			stS<<":";
			stmean<<":";
		}
	}
	bcal->Branch("adcN", b.adcN,sN.str().c_str());
	bcal->Branch("adcS", b.adcS,sS.str().c_str());
	bcal->Branch("tdcN", b.tdcN,stN.str().c_str());
	bcal->Branch("tdcS", b.tdcS,stS.str().c_str());
	bcal->Branch("tmean", b.tmean,stmean.str().c_str());
	bcal->Branch("adcsum",&b.adcsumN,"adcsumN/S:adcsumS/S");
	bcal->Branch("trig",&b.trig,"trig/I");
	bcal->Branch("Ntdctrig", &b.Ntdctrig,"Ntdctrig/I");
	bcal->Branch("tdctrig", b.tdctrig,"tdctrig[Ntdctrig]/F");
	bcal->Branch("tdctrigid", b.tdctrigid,"tdctrigid[Ntdctrig]/S");
	bcal->Branch("Nphotons", &b.Nphotons, "Nphotons/I");
	bcal->Branch("tphoton", b.tphoton, "tphoton[Nphotons]/F");
	bcal->Branch("Ephoton", b.Ephoton, "Ephoton[Nphotons]/F");
	bcal->Branch("tid", b.tid, "tid[Nphotons]/I");
	bcal->Branch("veto", &b.veto, "veto/F");
	bcal->Branch("tveto", &b.tveto, "tveto/F");

	beam_tree = new TTree("beam","Beam values from EPICS");
	beam_tree->Branch("beam",&beam.I2c24,"I2c24/F:x2c24:y2c24:I2c21:x2c21:y2c21:time/I");

	env_tree = new TTree("env","Environment values from EPICS");
	env_tree->Branch("env",&env.templ1,"templ1/F:templ2:templ3:humidityl1:humidityl2:humidityl3:time/I");

	// Unlock ROOT mutex
	pthread_mutex_unlock(&root_mutex);
	
	// Here we need to get the pedestals for this run so we have the
	// veto pedestal during event time. We get this from the 
	// DBCAL_ADCHit factory which reads them from the ROOT file.
	// Note that the pedestals are not actually read in yet since
	// we will get called before the brun routine of the DBCAL_ADCHit
	// factory. We don't really need the pedestals though, just
	// the pointer to the array where they will be when we do need them
	// during event time below.
	JFactory_base *jfac = eventLoop->GetFactory("DBCAL_ADCHit");
	if(jfac){
		DBCAL_ADCHit_factory *fac = (DBCAL_ADCHit_factory*)jfac;
		const float *gains1, *gains2; // unused
		float adc2gev; // unused
		fac->GetPedGains(ped_slot18, ped_slot19, gains1, gains2, adc2gev);
	}

	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t MyProcessor::evnt(JEventLoop *loop, int eventnumber)
{
	// If ROOT file is not open, then do nothing
	if(!file)return NOERROR;

	// If this is an EPICS event, we can add to the EPICS trees
	AddEPICSEvent(loop);

	// Get data
	vector<const DADC*> adcs;
	vector<const DTDC*> tdcs;
	vector<const DTrigger*> triggers;
	vector<const DBCAL_ADCHit*> adchits;
	vector<const DBCAL_TDCHit*> tdchits;
	vector<const DTAGGER_ETCoinc*> ethits;
	vector<const DTRIGGER_TDCHit*> trigtdcs;

	loop->Get(adcs);
	loop->Get(tdcs);
	loop->Get(triggers);
	loop->Get(adchits);
	loop->Get(tdchits);
	loop->Get(ethits);
	loop->Get(trigtdcs);

	// Initialize variables
	b.adcsumN = b.adcsumS = 0;
	b.trig=0;
	for(int i=0; i<18; i++)b.adcN[i] = b.adcS[i] = 0;
	for(int i=0; i<18; i++)b.tdcN[i] = b.tdcS[i] = -10000;

	// Trigger
	if(triggers.size()>0)b.trig=triggers[0]->latch;
	sort(trigtdcs.begin(), trigtdcs.end(), TrigTDCsort);
	b.Ntdctrig=0;
	float min_t = 1.0E6;
	int tdc_ref = 0; // holds TDC value of trigger that caused event
	int trig_ref = 1000000; // holds trigger bit used for ref (1-8)
	for(unsigned int i=0; i<trigtdcs.size(); i++){
		const DTRIGGER_TDCHit* tdc = trigtdcs[i];
		//if(fabs(tdc->t)>300.0)continue; // cut out obviously out-of-time hits
		if(b.Ntdctrig<MAX_TRIGTDCS){
			b.tdctrig[b.Ntdctrig] = tdc->t;
			b.tdctrigid[b.Ntdctrig++] = tdc->tid;
		}

		if(fabs(tdc->t)<min_t){
			min_t = fabs(tdc->t);
			tdc_ref=tdc->tdc;
			trig_ref = tdc->tid;
		}
		//if(b.Ntdctrig>=MAX_TRIGTDCS)break;
	}

	// BCAL ADC
	for(unsigned int i=0; i<adchits.size(); i++){
		const DBCAL_ADCHit *adchit = adchits[i];
		
		if(adchit->cid<1 || adchit->cid>18)continue;
		if(adchit->end==0){
			b.adcN[adchit->cid - 1] = adchit->adc-adchit->ped;
			b.adcsumN += adchit->adc-adchit->ped;
		}else{
			b.adcS[adchit->cid - 1] = adchit->adc-adchit->ped;
			b.adcsumS += adchit->adc-adchit->ped;
		}
	}
	
	// BCAL TDC
	for(unsigned int i=0; i<tdchits.size(); i++){
		const DBCAL_TDCHit *tdchit = tdchits[i];
		
		if(tdchit->cid<1 || tdchit->cid>18)continue;
		if(tdchit->end==0){
			b.tdcN[tdchit->cid - 1] = tdchit->t;
		}else{
			b.tdcS[tdchit->cid - 1] = tdchit->t;
		}
	}
	
	// BCAL mean times (timewalk corrected)
	for(int i=0; i<18; i++){
		double tn = b.tdcN[i];
		double ts = b.tdcS[i];
		double an = b.adcN[i];
		double as = b.adcS[i];
		
		// Apply timewalk corrections. We have functional forms from data
		// for cells 7,8, and 9. For the others, we use the cell 9 function
		// since there appears to be a rough correlation to the energy
		// deposition in the cell and cell 9 typically had less than 7 and 8.
		switch(i+1){
			case 7:
				tn = tn -(3.318633*exp(-0.007758*an)+-2.285192+(-0.000754*an));
				ts = ts -(2.209792*exp(-0.006442*as)+-2.561431+(-0.000607*as));
				break;
			case 8:
				tn = tn -(2.512104*exp(-0.004475*an)+-2.579980+(-0.000414*an));
				ts = ts -(2.396487*exp(-0.005468*as)+-2.408689+(-0.000530*as));
				break;
			case 9:
			default:
				tn = tn -(3.379855*exp(-0.008501*an)+-2.716681+(-0.001110*an));
				ts = ts -(3.284677*exp(-0.006599*an)+-2.976956+(-0.000848*an));
		}
		
		b.tmean[i] = (tn+ts)/2.0;
	}							

	// Tagger
	BCAL_TIME = b.tmean[8-1]; // used by TAGGERTDCsort
	sort(ethits.begin(), ethits.end(), TAGGERTDCsort);
	b.Nphotons = 0;
	for(unsigned int i=0; i<ethits.size(); i++){
		b.tphoton[b.Nphotons] = ethits[i]->time;
		b.Ephoton[b.Nphotons] = ethits[i]->energy;
		b.tid[b.Nphotons] = ethits[i]->T_id;
		b.Nphotons++;
		if(b.Nphotons>=MAX_PHOTONS)break;
	}
	
	// VETO
	b.veto = b.tveto = -1.0E6;
	for(unsigned int i=0; i<adcs.size(); i++){
		const DADC *adc = adcs[i];

		if(adc->crate != 0x1)continue;
		if(adc->slot != 19)continue;
		if(adc->channel != 29)continue;
		b.veto = adc->adc - ped_slot19[adc->channel];
	}
	for(unsigned int i=0; i<tdcs.size(); i++){
		const DTDC *tdc = tdcs[i];

		if(tdc->crate != 0x1)continue;
		if(tdc->slot != 11)continue;
		if(tdc->channel != 4)continue;
		b.tveto = (tdc->tdc - tdc_ref+100)*57.0E-12/1.0E-9; // convert to ns
	}

	// Fill bcal tree
	if(adchits.size()>0){
		// When we run with multiple threads, we should lock
		// a mutex so that only a single thread is accessing
		// the ROOT global section at a time. To make threading
		// efficient, we need to minimize the amount of time
		// spent with the mutex locked. Thus, all of the
		// loop->Get(...) calls should be done first and
		// pretty much just TTree Fill(...) methods
		// called while the mutex is locked.
		pthread_mutex_lock(&root_mutex);
		bcal->Fill();	
		pthread_mutex_unlock(&root_mutex);
	}

	return NOERROR;
}

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

//------------------
// fini
//------------------
jerror_t MyProcessor::fini(void)
{
	TThread::Lock();

	// Close file
	if(file){
		file->Write();
		file->Close();
		delete file;
	}
	file = NULL;

	TThread::UnLock();
	return NOERROR;
}

//------------------
// AddEPICSEvent
//------------------
void MyProcessor::AddEPICSEvent(JEventLoop *loop)
{
	vector<const DEPICS*> epics;
	loop->Get(epics);
	if(epics.size() == 0)return;

	// Initialize values
	beam.I2c24 = beam.x2c24 = beam.y2c24 = 0.0;
	beam.I2c21 = beam.x2c21 = beam.y2c21 = 0.0;
	env.templ1 = env.templ2 = env.templ3 = 0.0;
	env.humidityl1 = env.humidityl2 = env.humidityl3 = 0.0;
	
	
	// Loop over EPICS variables, copying the interesting ones
	for(unsigned int i=0; i<epics.size(); i++){
		const DEPICS *e = epics[i];
		beam.time = env.time = e->time;
		if(e->name == "hallb_IPM2C21A_CUR")beam.I2c21 = atof(e->value.c_str());
		if(e->name == "hallb_IPM2C21A_XPOS")beam.x2c21 = atof(e->value.c_str());
		if(e->name == "hallb_IPM2C21A_YPOS")beam.y2c21 = atof(e->value.c_str());
		if(e->name == "hallb_IPM2C24A_CUR")beam.I2c24 = atof(e->value.c_str());
		if(e->name == "hallb_IPM2C24A_XPOS")beam.x2c24 = atof(e->value.c_str());
		if(e->name == "hallb_IPM2C24A_YPOS")beam.y2c24 = atof(e->value.c_str());

		if(e->name == "hallb_l1_temp")env.templ1 = atof(e->value.c_str());
		if(e->name == "hallb_l2_temp")env.templ2 = atof(e->value.c_str());
		if(e->name == "hallb_l3_temp")env.templ3 = atof(e->value.c_str());
		if(e->name == "hallb_l1_hum")env.humidityl1 = atof(e->value.c_str());
		if(e->name == "hallb_l2_hum")env.humidityl2 = atof(e->value.c_str());
		if(e->name == "hallb_l3_hum")env.humidityl3 = atof(e->value.c_str());
	}
	// Check if there was beam current info
	if(beam.I2c21!=0.0 || beam.x2c21!=0.0 || beam.y2c21!=0.0){
		pthread_mutex_lock(&root_mutex);
		beam_tree->Fill();
		pthread_mutex_unlock(&root_mutex);
	}
	
	// Check if there was environment info
	if(env.templ1!=0.0 || env.templ2!=0.0 || env.templ3!=0.0){
		pthread_mutex_lock(&root_mutex);
		env_tree->Fill();
		pthread_mutex_unlock(&root_mutex);
	}

}