// $Id$
//
//    File: DBCALTimeSpectrum_factory.cc
// Created: Mon Aug  1 08:55:45 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 <TRACKING/DMCTrajectoryPoint.h>

#include "DBCALTimeSpectrum_factory.h"
using namespace jana;

// These are auto-generated from timewalk_calibrate_max.C for
// the specific summing scheme implemented here.
static double BCAL_TimeWalkUp(double fADC, int layer);
static double BCAL_TimeWalkDn(double fADC, int layer);
#include "timewalk_FINE.h"

#define M_2PI 6.283185307

// see slides shown at 7/21/2011 BCAL segmentation meeting
// n.b. the number of photons/MeV/side for the actual fibers
// has changed from 75 to 290/2=145 (long range part).
// To use the old 75 number, set the MEV_PER_PE configuration
// parameter to 0.668.
double MeV_per_PE = 0.346;


//------------------
// init
//------------------
jerror_t DBCALTimeSpectrum_factory::init(void)
{
	// Configuration parameters
	gPARMS->SetDefaultParameter("MEV_PER_PE", MeV_per_PE, "MeV per photo electron produced in fiber and seen by SiPM.");

	// We want to create a DBCALTimeSpectrum object for each
	// SiPM cell (layer/sector). We keep these around so that
	// we don't have to keep reallocating memory for them
	// (and the DHistogram objects they contain!)
	
	// Inner
	for(int layer=1; layer<=6; layer++){
		for(int sector=1; sector<=4; sector++){
			pair<int,int> idx = pair<int,int>(layer, sector);
			bcaltimespectra[idx] = new DBCALTimeSpectrum(layer, sector);
			vbcaltimespectra.push_back(bcaltimespectra[idx]);
		}
	}

	// Outer
	for(int layer=7; layer<=10; layer++){
		for(int sector=1; sector<=4; sector++){
			pair<int,int> idx = pair<int,int>(layer, sector);
			bcaltimespectra[idx] = new DBCALTimeSpectrum(layer, sector);
			vbcaltimespectra.push_back(bcaltimespectra[idx]);
		}
	}

	// By setting the NOT_OBJECT_OWNER flag, the JANA framework
	// will not delete the objects in _data at the start of the
	// event allowing us to keep them around. (n.b. using the
	// PERSISTANT flag would cause JANA to skip calling out evnt
	// method so we don't use that!)
	flags = NOT_OBJECT_OWNER;

	return NOERROR;
}

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

//------------------
// evnt
//------------------
jerror_t DBCALTimeSpectrum_factory::evnt(JEventLoop *loop, int eventnumber)
{
	// Re-copy pointers into _data (see note on NOT_OBJECT_OWNER in 
	// init above)
	_data = vbcaltimespectra;

	// Clear all DBCALTimeSpectrum objects
	for(unsigned int i=0; i<_data.size(); i++)_data[i]->Reset();

	// Get all DMCTrajectoryPoint objects for this event
	vector<const DMCTrajectoryPoint*> mctrajpoints;
	loop->Get(mctrajpoints);

	double ATTEN_LENGTH = 300.0;
	double C_EFFECTIVE  = 16.75;
		
	// Loop over trajectory point objects and use them to fill the
	// DBCALTimeSpectrum attenuated energy histograms.
	double Emax = 0.0;
	double theta_thrown = 0.0;
	for(unsigned int i=0; i<mctrajpoints.size(); i++){
		const DMCTrajectoryPoint *traj = mctrajpoints[i];

		// More convenient to work with single letter variables
		double X = traj->x;
		double Y = traj->y;
		double Z = traj->z;
		double my_dE = traj->dE;

		// n.b.! The step stores dE in GeV, but the toStrings method
		// converts it to MeV for printing. Thus, janaroot already has
		// the value in MeV, but when we deal with the DMCTrajectoryPoint
		// objects in C++, they are in GeV and must be converted.
		my_dE *= 1000.0;

		// record largest energy particle's energy
		// and assume that is incident energy.
		if(traj->E>Emax && traj->part<=3){
			Emax = traj->E;
			double p = sqrt(traj->px*traj->px + traj->py*traj->py + traj->pz*traj->pz);
			theta_thrown = acos(traj->pz/p);
		}
		
		// Make sure step is inside BCAL sensitive volume
		if(Z<17.0 || Z>407.0)continue;
		
		double R = sqrt(X*X + Y*Y);
		if(R<64.3 || R>86.172)continue;
		
		// Find SiPM layer
		int layer=1;
		if(R>66.300) layer++;
		if(R>68.300) layer++;
		if(R>70.300) layer++;
		if(R>72.300) layer++;
		if(R>74.300) layer++;
		if(R>76.300) layer++;
		if(R>78.768) layer++;
		if(R>81.236) layer++;
		if(R>83.704) layer++;

		double phi = atan2(Y,X);
		if(phi<0.0)phi += M_2PI;
		
		double phi_one_module = M_2PI/48.0;
		double module = 1.0 + floor(phi/phi_one_module);
		double phi_rel = phi - (module-1.0)*phi_one_module;
		double phi_one_cell = phi_one_module/4.0;
		int sector = 1 + (int)floor(phi_rel/phi_one_cell);
		
		// Get map index based on layer and sector
		pair<int,int> idx = pair<int,int>(layer,sector);
		
		if(bcaltimespectra.find(idx) == bcaltimespectra.end()){
			cout<<"Not found: module:"<<module<<" layer:"<<layer<<" sector:"<<sector<<endl;
		}else{
			double dist_up = Z - 17.0;
			double dist_dn = 407.0 - Z;
			
			double Eup = my_dE*exp(-dist_up/ATTEN_LENGTH);
			double Edn = my_dE*exp(-dist_dn/ATTEN_LENGTH);
			double t_up = traj->t*1.0E9 + dist_up/C_EFFECTIVE;
			double t_dn = traj->t*1.0E9 + dist_dn/C_EFFECTIVE;
			
			DBCALTimeSpectrum *spectrum = (bcaltimespectra.find(idx))->second;
			spectrum->Eup.Fill(t_up, Eup);
			spectrum->Edn.Fill(t_dn, Edn);
			spectrum->Etot += my_dE;

//if(traj->dE!=0.0)_DBG_<<"tup:"<<t_up<<" tdn:"<<t_dn<<" Eup:"<<Eup<<" Edn:"<<Edn<<" layer:"<<layer<<" sector:"<<sector<<" module:"<<module<<endl;
		}
	}
	
	// Loop over all cells
	for(unsigned int i=0; i<_data.size(); i++){
		
		DBCALTimeSpectrum *spectrum = _data[i];
	
		// Smear each of the attenuated energy spectra due to sampling
		// fluctuations.
		spectrum->SamplingSmear(Emax, theta_thrown);

		// Include Poisson statistics
		spectrum->PoissonSmear();
		
		// Include time jitter of detector
		spectrum->DetectorJitterSmear();

		// Add Dark hits to each SiPM
		spectrum->AddDarkPulses();

		// Convolute the smeared, attenuated energy spectrum
		// with the electronic pulse shape
		spectrum->Convolute();
		
		// Apply threshold.
		double thresh_mV = 44.7;
		spectrum->FindTimes(thresh_mV);
	}

	return NOERROR;
}


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

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