// $Id$
//
//    File: DBCALTimeSpectrum_factory_2006BeamTest.cc
// Created: Thu Aug  30 9:34:05 EDT 2011
// Creator: davidl (on Linux ifarm1102 2.6.18-128.7.1.el5 x86_64)
//


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

#include "DBCALTimeSpectrum_factory.h"
#include "DBCALTimeSpectrum_factory_2006BeamTest.h"
using namespace jana;

// These are copied from a version of timewalk.cc which is
// auto-generated from timewalk_calibrate.C for the specific
// summing scheme implemented here.
static double BCAL_TimeWalkUp(double fADC, int layer);
static double BCAL_TimeWalkDn(double fADC, int layer);

//------------------
// init
//------------------
jerror_t DBCALTimeSpectrum_factory_2006BeamTest::init(void)
{
	/// Fill in the sum_map. This is used later to produce summed
	/// cells.

	// It's easiest just to loop over all SiPM layer/sector combos
	// and index the map with the id causing it to be created if it
	// doesn't already exist.
	for(int layer=1; layer<=10; layer++){ 
		for(int sector=1; sector<=4; sector++){
		
			int idx = GetSummedCell(layer, sector);

			sum_map[idx]; // This is all that's need to ensure an entry in sum_map exists
		}
	}

	return NOERROR;
}


//------------------
// GetSummedLayer
//------------------
int DBCALTimeSpectrum_factory_2006BeamTest::GetSummedLayer(int layer, int sector)
{
	switch(layer){
		case 1:
		case 2:
			return 1;
		case 3:
		case 4:
			return 2;
		case 5:
		case 6:
			return 3;
		case 7:
			return 4;
		case 8:
			if(sector==2 || sector==3)return 5;
			return 4;
		case 9:
			if(sector==2 || sector==3)return 6;
			return 6;
		case 10:
			if(sector==2 || sector==3)return 7;
			return 6;
	}

	return 0;
}

//------------------
// GetSummedSector
//------------------
int DBCALTimeSpectrum_factory_2006BeamTest::GetSummedSector(int layer, int sector)
{
	switch(sector){
		case 1:
			return 1;
		case 2:
		case 3:
			return 2;
		case 4:
			return 3;
	}
	
	return 0;
}

//------------------
// GetSummedCell
//------------------
int DBCALTimeSpectrum_factory_2006BeamTest::GetSummedCell(int layer, int sector)
{
	// Return the index of the summed cell given the layer
	// and sector. 
	
	int summed_layer = GetSummedLayer(layer, sector);
	int summed_sector = GetSummedSector(layer, sector);
	
	return summed_layer*10 + summed_sector;
}

//------------------
// GetSummedLayerAndSectorFromCell
//------------------
void DBCALTimeSpectrum_factory_2006BeamTest::GetSummedLayerAndSectorFromCell(int cell, int &summed_layer, int &summed_sector)
{
	// Return the summed layer and summed sector based on the cell number
	
	summed_layer = cell/10;
	summed_sector = cell%10;
}

//------------------
// brun
//------------------
jerror_t DBCALTimeSpectrum_factory_2006BeamTest::brun(jana::JEventLoop *eventLoop, int runnumber)
{
	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t DBCALTimeSpectrum_factory_2006BeamTest::evnt(JEventLoop *loop, int eventnumber)
{
	// Get the single SiPM spectra
	vector<const DBCALTimeSpectrum *> spectra;
	
	// Jump through hoops to get to untagged factory
	DBCALTimeSpectrum_factory *fac = (DBCALTimeSpectrum_factory*)loop->GetFactory("DBCALTimeSpectrum","",false);
	fac->Get(spectra);

	// Clear all sum vectors
	map<int, vector<const DBCALTimeSpectrum*> >::iterator iter;
	for(iter=sum_map.begin(); iter!=sum_map.end(); iter++){
		iter->second.clear();
	}

	// Loop over single SiPM spectra and sort them into 
	// lists according to the summing scheme
	for(unsigned int i=0; i<spectra.size(); i++){
		const DBCALTimeSpectrum *spectrum = spectra[i];
		
		int idx = GetSummedCell(spectrum->layer, spectrum->sector);
		sum_map[idx].push_back(spectrum);
	}

	// Loop over all summed cells and create DBCALTimeSpectrum objects
	// for each
	for(iter=sum_map.begin(); iter!=sum_map.end(); iter++){
		
		int summed_layer, summed_sector;
		GetSummedLayerAndSectorFromCell(iter->first, summed_layer, summed_sector);
		
		double thresh_mV = 44.7;
		DBCALTimeSpectrum *spectrum = new DBCALTimeSpectrum(thresh_mV, iter->second, summed_layer, summed_sector);
		
		// Overwrite corrected times so we can apply timewalk
		// corrections for the summed cells
		spectrum->tup_corrected = spectrum->tup;
		if(spectrum->tup>-1000.0){
			spectrum->tup_corrected -= BCAL_TimeWalkUp(spectrum->fADC_up, summed_layer);
		}
		spectrum->tdn_corrected = spectrum->tdn;
		if(spectrum->tdn>-1000.0){
			spectrum->tdn_corrected -= BCAL_TimeWalkDn(spectrum->fADC_dn, summed_layer);
		}
		
		_data.push_back(spectrum);
	}

	return NOERROR;
}

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

//------------------
// fini
//------------------
jerror_t DBCALTimeSpectrum_factory_2006BeamTest::fini(void)
{
	return NOERROR;
}


#define NLAYERS_BCAL 10
double BCAL_TimeWalkUp(double fADC, int layer)
{
	double tw_par[NLAYERS_BCAL][4] ={
		{18.9481, 0.394631, 0.0747397, -1.26256},	// layer 1
		{18.8913, 0.0792789, 0.0171958, -1.05491},	// layer 2
		{18.5248, 0.666298, 0.0718911, -1.20172},	// layer 3
		{19.3237, 0.38947, 0.0734877, -1.2575},	// layer 4
		{18.9397, 0.260441, 0.0366492, -1.10454},	// layer 5
		{19.3164, 0.0756342, 0.0163968, -1.05227},	// layer 6
		{19.8672, 0.0983631, 0.0266604, -1.08922},	// layer 7
		{40, 0.8, 0.064, -1.2},	// layer 8
		{40, 0.8, 0.064, -1.2},	// layer 9
		{40, 0.8, 0.064, -1.2}	// layer 10
		};

	layer = layer<1 ? 1:layer>NLAYERS_BCAL ? NLAYERS_BCAL:layer;
	double *p = tw_par[layer-1];
	return p[0]+p[1]/(pow(fADC,p[2])+p[3]);
}


double BCAL_TimeWalkDn(double fADC, int layer)
{
	double tw_par[NLAYERS_BCAL][4] ={
		{18.8394, 0.0170437, 0.00413264, -1.01302},	// layer 1
		{18.1125, 0.04527, 0.00560642, -1.01415},	// layer 2
		{18.3564, 0.0487526, 0.00652681, -1.0172},	// layer 3
		{19.1427, 0.0199204, 0.00452909, -1.0142},	// layer 4
		{19.8346, 0.0338277, 0.0104048, -1.03418},	// layer 5
		{18.6599, 0.0481108, 0.00645093, -1.01704},	// layer 6
		{20.2037, 0.0935863, 0.0301207, -1.10329},	// layer 7
		{17, 0.2, 0.034, -1.1},	// layer 8
		{17, 0.2, 0.034, -1.1},	// layer 9
		{17, 0.2, 0.034, -1.1}	// layer 10
		};

	layer = layer<1 ? 1:layer>NLAYERS_BCAL ? NLAYERS_BCAL:layer;
	double *p = tw_par[layer-1];
	return p[0]+p[1]/(pow(fADC,p[2])+p[3]);
}