// $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/DMCThrown.h>
#include <TRACKING/DMCTrajectoryPoint.h>

#include "DBCALTimeSpectrum_factory.h"
using namespace jana;

// These are auto-generated from timewalk_calibrate.C 
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;

// Threshold for TDC in mV
double THRESH_TDC_mV = 44.7;

// This value can be used to apply additional gain
// to the signal. It does not correspond directly
// to a physical amplifier, but will be used to
// scale both the fADC and TDC signals
double ADDITIONAL_GAIN=1.0;

// Gain of the second stage pre-amp used for TDC only
// The pulse shape used here represents an amplification factor
// of 10 for the TDC leg. Elton claims this is way too high
// and will be lowered to probably no more than 5. 
double PREAMP_GAIN_TDC=5.0;

// Sampling fluctuation parameters. Normally, these
// come from the calibDB. The floor term BCAL_SAMPLINGCOEFB
// is set to 0.013 in sim-recon, but it's unclear
// what the value should be. Also, the floor term is
// added in quadrature here while it is added linearly
// in sim-recon
double BCAL_SAMPLINGCOEFA = 0.042;
double BCAL_SAMPLINGCOEFB = 0.0; // See slide 14 from June 17th BCAL segmentation meeting slides

// Flags used to turn off various effects
bool APPLY_SAMPLING_SMEAR = true;
bool APPLY_PHOTOSTATISTICS = true;
bool APPLY_TIME_JITTER = true;
bool APPLY_DARK_HITS = true;
bool APPLY_ELECTRONIC_NOISE = true;

//------------------
// 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.");
	gPARMS->SetDefaultParameter("THRESH_TDC_mV", THRESH_TDC_mV, "Threshold applied to electronic pulse to determine TDC time.");
	gPARMS->SetDefaultParameter("ADDITIONAL_GAIN", ADDITIONAL_GAIN, "Additional gain applied to both fADC and TDC");
	gPARMS->SetDefaultParameter("PREAMP_GAIN_TDC", PREAMP_GAIN_TDC, "Second stage gain for just TDC");

	gPARMS->SetDefaultParameter("BCAL_SAMPLINGCOEFA", BCAL_SAMPLINGCOEFA, "Sampling fluctuations stochastic term");
	gPARMS->SetDefaultParameter("BCAL_SAMPLINGCOEFB", BCAL_SAMPLINGCOEFB, "Sampling fluctuations floor term");

	gPARMS->SetDefaultParameter("APPLY_SAMPLING_SMEAR", APPLY_SAMPLING_SMEAR, "Switch to turn on/off sampling fluctuations");
	gPARMS->SetDefaultParameter("APPLY_PHOTOSTATISTICS", APPLY_PHOTOSTATISTICS, "Switch to turn on/off photostatistics");
	gPARMS->SetDefaultParameter("APPLY_TIME_JITTER", APPLY_TIME_JITTER, "Switch to turn on/off time jitter in SiPM");
	gPARMS->SetDefaultParameter("APPLY_DARK_HITS", APPLY_DARK_HITS, "Switch to turn on/off dark hit generation");
	gPARMS->SetDefaultParameter("APPLY_ELECTRONIC_NOISE", APPLY_ELECTRONIC_NOISE, "Switch to turn on/off electronic noise");

	// 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)
{
	const DMCThrown *thrown;
	loop->GetSingle(thrown);

	// 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 = thrown->energy();
	double theta_thrown = thrown->momentum().Theta();
	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;
		}
	}
	
	// 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.
		if(APPLY_SAMPLING_SMEAR)spectrum->SamplingSmear(Emax, theta_thrown);

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

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

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

	return NOERROR;
}


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

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