// $Id$
//
// Author: David Lawrence  August 27, 2006
//
//
// DEventSourceROOT_mc methods
//

#include <stdlib.h>

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

#include <TThread.h>

#include <JANA/JApplication.h>
#include <JANA/JFactory_base.h>
#include <JANA/JFactory.h>
#include <JANA/JEventLoop.h>
#include <JANA/JEvent.h>

#include "DEventSourceROOT_mc.h"
#include "BCAL/DBCAL_ADCHit_factory.h"

// Sorting functions used to put TDC and ADC hits in geographic order
bool ADCsort(DADC* const &thit1, DADC* const &thit2) {
	if(thit1->crate != thit2->crate)return thit1->crate < thit2->crate;
	if(thit1->slot != thit2->slot)return thit1->slot < thit2->slot;
	return thit1->channel < thit2->channel;
}
bool TDCsort(DTDC* const &thit1, DTDC* const &thit2) {
	if(thit1->crate != thit2->crate)return thit1->crate < thit2->crate;
	if(thit1->slot != thit2->slot)return thit1->slot < thit2->slot;
	return thit1->channel < thit2->channel;
}


//----------------
// Constructor
//----------------
DEventSourceROOT_mc::DEventSourceROOT_mc(const char* source_name):JEventSource(source_name)
{
	// Open ROOT file and get pointer to TTree named h999
	TThread::Lock();
	tree = NULL;
	file = new TFile(source_name);
	if(!file)return;
	
	// Get the TTree
	TObject *obj = (TObject*)file->Get("h999");
	if(obj){
		tree = dynamic_cast<TTree*>(obj);
	}
	
	// Set up branch addresses
	if(tree){
		tree->SetBranchAddress("Nhits", &Nhits);
		tree->SetBranchAddress("Zin", &Zin);
		tree->SetBranchAddress("Nh", &Nh);
		tree->SetBranchAddress("Nv", &Nv);
		tree->SetBranchAddress("Escih", &Escih);
	}
	
	
	TThread::UnLock();
}

//----------------
// Destructor
//----------------
DEventSourceROOT_mc::~DEventSourceROOT_mc()
{
	// Close file
	if(tree) delete tree;
	if(file) delete file;
}

//----------------
// GetEvent
//----------------
jerror_t DEventSourceROOT_mc::GetEvent(JEvent &event)
{
	/// Implementation of JEventSource::GetEvent function

	// For this source we do not do optimize too much. Each
	// request will read in the event from the file and
	// grab out the info. Specifically, if TDC and ADC info
	// are requested, it will read the event in twice from
	// the ROOT file. While this wastes some time, it is simple
	// to write and this source is not expected to be used too
	// terribly much.
	
	// Event number is not kept in file. Just keep a counter
	int event_number = ++Nevents_read;
	int run_number = 100;
	
	// Check that event is in range
	TThread::Lock();
	int Nentries = tree->GetEntries();
	TThread::UnLock();
	if(Nevents_read>Nentries)return NO_MORE_EVENTS_IN_SOURCE;
	
	// Copy the reference info into the JEvent object
	event.SetJEventSource(this);
	event.SetEventNumber(event_number);
	event.SetRunNumber(run_number);
	event.SetRef((void*)(Nevents_read-1));

	return NOERROR;
}

//----------------
// FreeEvent
//----------------
void DEventSourceROOT_mc::FreeEvent(JEvent &event)
{
	// Don't need to free anything
}

//----------------
// GetObjects
//----------------
jerror_t DEventSourceROOT_mc::GetObjects(JEvent &event, JFactory_base *factory)
{
	/// This gets called through the virtual method of the
	/// JEventSource base class. It creates the objects of the type
	/// on which factory is based.
	
	// As stated above, this is not done in the most efficient way.
	// Read in the entire event and extract the desired info from it.

	// We must have a factory to hold the data
	if(!factory)throw RESOURCE_UNAVAILABLE;
	
	// Don't support tagged factories
	if(strcmp(factory->Tag(), ""))return OBJECT_NOT_AVAILABLE;
	
	// The ref field of the JEvent is just the event number.(starting from 0)
	int eventNo = (int)event.GetRef();
	if(!tree)throw RESOURCE_UNAVAILABLE;
	
	// Figure out what type of data is being requested
	JFactory<DADC> *fac_adc = dynamic_cast<JFactory<DADC>*>(factory);
	if(fac_adc)return Extract_DADC(eventNo, fac_adc, event.GetJEventLoop());

	JFactory<DTDC> *fac_tdc = dynamic_cast<JFactory<DTDC>*>(factory);
	if(fac_tdc)return Extract_DTDC(eventNo, fac_tdc);

	JFactory<DTrigger> *fac_trig = dynamic_cast<JFactory<DTrigger>*>(factory);
	if(fac_trig)return Extract_DTrigger(eventNo, fac_trig);
		
	return OBJECT_NOT_AVAILABLE;
}



//----------------
// Extract_DTrigger
//----------------
jerror_t DEventSourceROOT_mc::Extract_DTrigger(int eventNo, JFactory<DTrigger> *fac)
{
	/// Get trigger information from Trigger Supervisor bank

	if(!fac)return OBJECT_NOT_AVAILABLE;
	vector<DTrigger*> data;
	
	DTrigger *trig = new DTrigger;
	trig->latch = 0x10;
	trig->live1 = 10*eventNo;
	trig->live2 = 10*eventNo;
	data.push_back(trig);
	
	fac->CopyTo(data);

	return NOERROR;
}


//----------------
// Extract_DADC
//----------------
jerror_t DEventSourceROOT_mc::Extract_DADC(int eventNo, JFactory<DADC> *fac, JEventLoop *loop)
{
	if(!fac)return OBJECT_NOT_AVAILABLE;
	vector<DADC*> data;

	// Read in the event
	TThread::Lock();
	tree->GetEntry(eventNo);
	
	// Get conversion values from DBCAL_ADCHit factory so we can convert
	// energy to ADC in such a way that it will be properly converted back
	DBCAL_ADCHit_factory *adchit_fac = NULL;
	JFactory_base *fac_base = loop->GetFactory("DBCAL_ADCHit");
	if(fac_base)adchit_fac = dynamic_cast<DBCAL_ADCHit_factory*>(fac_base);
	const int *ped_slot18;
	const int *ped_slot19;
	const float *gain_slot18;
	const float *gain_slot19;
	float ADC_TO_GEV;
	if(adchit_fac)adchit_fac->GetPedGains(ped_slot18, ped_slot19, gain_slot18, gain_slot19, ADC_TO_GEV);
	
	// Loop over sections hit, creating DADC hits for both ends of the BCAL
	float z = Zin[0] + (390.0/2.0);
	for(int i=0; i<(int)Nhits[0]; i++){
		// Looking at the beam right end, the section ids are:
		// 
		//     -------------------------------
		//     |  1 |  2 |  3 |  4 |  5 |  6 | 
		//     -------------------------------
		//     |  7 |  8 |  9 | 10 | 11 | 12 | 
		//     -------------------------------
		//     | 13 | 14 | 15 | 16 | 17 | 18 | 
		//     -------------------------------
		//
		// ADCs are in slots 18(BCALN) and 19(BCALS).
		// The first 18 channels of each (channels 0-17)
		// correspond to the id in the above chart.
		
		int id = Nv[i] + (Nh[i]-1)*6;
		int channel = id-1;
		float E = Escih[i];
		
		// Attenuate to ends
		float lambda = 300.0; // effective attenuation length of BCAL is 300cm
		float L = 390.0; // length of modules is 390cm
		float E_left = E*exp(-z/lambda);
		float E_right = E*exp(-(L-z)/lambda);
		
		// Convert to ADC values
		int adc_left = (int)(E_left*20000.0) + 100;
		int adc_right = (int)(E_right*20000.0) + 100;
		if(adchit_fac){
			// use DBCAL_ADCHit conversion values to convert to ADC units
			adc_left = (int)(E_left/gain_slot19[channel]/ADC_TO_GEV) + ped_slot19[channel];
			adc_right = (int)(E_right/gain_slot18[channel]/ADC_TO_GEV) + ped_slot18[channel];
		}
		
		// BCALN (beam right)
		DADC *adc = new DADC;
		adc->crate = 0x1;
		adc->slot = 18;
		adc->channel = channel;
		adc->adc = adc_right;
		data.push_back(adc);

		// BCALS (beam left)
		adc = new DADC;
		adc->crate = 0x1;
		adc->slot = 19;
		adc->channel = channel;
		adc->adc = adc_left;
		data.push_back(adc);
	}
	
	TThread::UnLock();

	sort(data.begin(), data.end(), ADCsort);
	fac->CopyTo(data);

	return NOERROR;
}

//----------------
// Extract_DTDC
//----------------
jerror_t DEventSourceROOT_mc::Extract_DTDC(int eventNo, JFactory<DTDC> *fac)
{
	if(!fac)return OBJECT_NOT_AVAILABLE;
	vector<DTDC*> data;
	
	// Read in the event
	TThread::Lock();
	tree->GetEntry(eventNo);
	
	// The F1TDC does not automatically subtract the TDC stop time. A
	// channel must be used for reference. We will use the trigger
	// signals since they should be in time with the TDC stop. Here,
	// generate a random value for the trigger time.
	DTDC *trig_tdc = new DTDC;
	trig_tdc->crate = 0x1;
	trig_tdc->slot = 11;
	trig_tdc->channel = 20; // 5th bit of second connector. channels start at 0
	trig_tdc->tdc = random()&0x3ff;
	data.push_back(trig_tdc);
	
	// Loop over sections hit, creating DTDC hits for both ends of the BCAL
	float z = Zin[0] + (390.0/2.0);
	for(int i=0; i<(int)Nhits[0]; i++){
		// See Extract_DADC for description of numbering scheme.
		//
		// TDCs are in slots 10 and 11. Because of how we need to
		// group discriminator channels for the trigger, the TDC
		// channel numbers do not line up with the ADC ones.
		// The first two connectors in the TDC in slot 10 will
		// contain BCAL TDC channels for all sections except those
		// with id 1 and 13. Those sections will go into the
		// first 4 channels (2 from BCALN and 2 from BCALS) of the
		// TDC in slot 11.
		//
		// To add some realism, we impose a threshold here on the
		// attenuated energy being over 1MeV at the end. Note
		// that there is about a 60% of the single from the middle
		// of the BCAL to the ends so this means an effective threshold
		// of 17MeV for showers starting in the middle.
		
		int id = Nv[i] + (Nh[i]-1)*6;
		int slot = 10;
		int channel = id-1;
		float E = Escih[i];
		
		// Attenuate to ends
		float lambda = 300.0; // effective attenuation length of BCAL is 300cm
		float L = 390.0; // length of modules is 390cm
		float E_left = E*exp(-z/lambda);
		float E_right = E*exp(-(L-z)/lambda);
		
		// BCALN (beam right)
		if(E_right>0.001){
			// assume effective speed of light in BCAL is 0.9c
			// TDC resolution is about 60ps
			float c = 2.998E10;
			float t = (L-z)/(0.9*c); // in seconds
			float t_c = t/60.0E-12;
		
			switch(id){
				case  1:	channel = 0;	slot = 11;	break;
				case 13:	channel = 1;	slot = 11;	break;
			}
		
			DTDC *tdc = new DTDC;
			tdc->crate = 0x1;
			tdc->slot = slot;
			tdc->channel = channel;
			tdc->tdc = (int)t_c + 400 + trig_tdc->tdc + (random()&0x7); // offset=400
			data.push_back(tdc);
		}

		// BCALS (beam left)
		if(E_left>0.001){
			// assume effective speed of light in BCAL is 0.9c
			// TDC resolution is about 60ps
			float c = 2.998E10;
			float t = z/(0.9*c); // in seconds
			float t_c = t/60.0E-12;
		
			switch(id){
				case  1:	channel = 2;	slot = 11;	break;
				case 13:	channel = 3;	slot = 11;	break;
				default: channel += 16;
			}
		
			DTDC *tdc = new DTDC;
			tdc->crate = 0x1;
			tdc->slot = slot;
			tdc->channel = channel;
			tdc->tdc = (int)t_c + 400 + trig_tdc->tdc + (random()&0x7); // offset=400
			data.push_back(tdc);
		}
	}
	TThread::UnLock();

	sort(data.begin(), data.end(), TDCsort);
	fac->CopyTo(data);

	return NOERROR;
}