// $Id$
//
//    File: DBCALTimeSpectrum.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 <cmath>
using namespace std;

#include "DBCALTimeSpectrum.h"

// These are defined in timewalk.cc which is
// auto-generated from timewalk_calibrate.C
double BCAL_TimeWalkUp(double fADC, int layer);
double BCAL_TimeWalkDn(double fADC, int layer);


//---------------
// DBCALTimeSpectrum
//---------------
DBCALTimeSpectrum::DBCALTimeSpectrum(int layer, int sector):
			layer(layer),
			sector(sector),
			Eup(NBINS_E, E_LO, E_HI),
			Edn(NBINS_E, E_LO, E_HI),
			electronic_up(NBINS_ELECTRONIC, ELECTRONIC_LO, ELECTRONIC_HI),
			electronic_dn(NBINS_ELECTRONIC, ELECTRONIC_LO, ELECTRONIC_HI),
			rnd(0)
{
	Reset();
}

//---------------
// DBCALTimeSpectrum
//---------------
DBCALTimeSpectrum::DBCALTimeSpectrum(double thresh_mV, vector<const DBCALTimeSpectrum*> &spectra, int layer, int sector):
			layer(layer),
			sector(sector),
			Eup(NBINS_E, E_LO, E_HI),
			Edn(NBINS_E, E_LO, E_HI),
			electronic_up(NBINS_ELECTRONIC, ELECTRONIC_LO, ELECTRONIC_HI),
			electronic_dn(NBINS_ELECTRONIC, ELECTRONIC_LO, ELECTRONIC_HI),
			rnd(0)
{
	/// This constructor is used to sum together multiple SiPM
	/// timing distributions into a single distribution. It is
	/// assumed that the inputs already have filled histograms.
	/// The FindTimes() method is called automatically after
	/// the summing is complete.

	Reset();
	
	for(unsigned int i=0; i<spectra.size(); i++){
		const DBCALTimeSpectrum *spectrum = spectra[i];
		
		Eup.Add(&spectrum->Eup);
		Edn.Add(&spectrum->Edn);
		electronic_up.Add(&spectrum->electronic_up);
		electronic_dn.Add(&spectrum->electronic_dn);
		
		Etot += spectrum->Etot;
		fADC_up += spectrum->fADC_up;
		fADC_dn += spectrum->fADC_dn;
	}

	geometric_mean = sqrt((double)fADC_up*(double)fADC_dn);
	FindTimes(thresh_mV);
}

//---------------
// AddDarkPulses
//---------------
void DBCALTimeSpectrum::AddDarkPulses(void)
{
	// Add random dark hits to every event in the given histogram.

	double tmin = ELECTRONIC_LO;
	double tmax = ELECTRONIC_HI;

	double BCAL_DARKRATE_GHZ = 0.0176;
	double BCAL_XTALK_FRACT = 0.157;
	double BCAL_INTWINDOW_NS = tmax - tmin;
	double MeV_per_PE = 0.668; // see slides shown at 7/21/2011 BCAL segmentation meeting

	// Get dark pulses plus cross talk
	int darkPulse_up = (int)floor(0.5+rnd.Poisson( BCAL_DARKRATE_GHZ* BCAL_INTWINDOW_NS ));
	int darkPulse_dn = (int)floor(0.5+rnd.Poisson( BCAL_DARKRATE_GHZ* BCAL_INTWINDOW_NS ));
	int xTalk_up = (int)floor(0.5+rnd.Poisson( (double)darkPulse_up * BCAL_XTALK_FRACT ));
	int xTalk_dn = (int)floor(0.5+rnd.Poisson( (double)darkPulse_dn * BCAL_XTALK_FRACT ));

	// Convert from PE to MeV
	double dE_up = (double)( xTalk_up + darkPulse_up )*MeV_per_PE;
	double dE_dn = (double)( xTalk_dn + darkPulse_dn )*MeV_per_PE;

	// Randomly place the energy in time
	if(dE_up>0.0){
		double s  = rnd.Rndm();
		double t  = tmin + s*BCAL_INTWINDOW_NS;
		Eup.Fill(t, dE_up);
	}

	if(dE_dn>0.0){
		double s  = rnd.Rndm();
		double t  = tmin + s*BCAL_INTWINDOW_NS;
		Edn.Fill(t, dE_dn);
	}
	
}

//---------------
// Convolute
//---------------
void DBCALTimeSpectrum::Convolute(void)
{
	// The following will initialize a table that will be used
	// for convoluting the electronic pulse_shape function. It
	// is done this way to speed things up below. It is done here
	// since the table that is calculated is NBINS_E x NBINS_ELECTRONIC
	// in format. The possibility of multiple threads adds a little
	// more complication.
	static double *pulse_shape_pre=NULL;
	if(pulse_shape_pre==NULL){
	
		static pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
		pthread_mutex_lock(&mutex);

		if(pulse_shape_pre==NULL){
			cout<<"Pre-evaluating pulse_shape function ..."<<endl;
			pulse_shape_pre = new double[NBINS_E*NBINS_ELECTRONIC];
			
			// Get electronic pulse shape
			TSpline *pulse_shape = MakeTSpline();

			// Conversion factor between attenuated MeV and mV
			double mV_per_MeV = 8.82; // slides from BCAL segmentation meeting 7/21/2011 
	
			// Normalize pulse shape between t=15ns and t=100ns, excluding the
			// pre and post bumps. We scale by the amplitude to give the pulse
			// shape a height of 1 so that when multiplied by mV_per_MeV, the
			// amplitude will again be in mV.
			//double norm = pulse_shape->GetMinimum(15.0, 100.0);
			double norm = -1020.80; // mV not trivial to find minimum with TSpline3

			// time shift pulse shape just to bring it in frame better for zoomed in picture
			double toffset = 6.0;

			for(int ibin=1; ibin<=NBINS_E; ibin++){

				double t0 = Eup.GetBinCenter(ibin);

				for(int jbin=1; jbin<=NBINS_ELECTRONIC; jbin++){
					double t = electronic_up.GetBinCenter(jbin);
					int idx = jbin + (ibin-1)*NBINS_ELECTRONIC;
					pulse_shape_pre[idx] = (mV_per_MeV/norm)*pulse_shape->Eval(t-t0+toffset);
				}
			}
			delete pulse_shape;
			cout<<"Done."<<endl;
		}
		pthread_mutex_unlock(&mutex);
	}

	// Convolute timing histos with pulse shape
	for(int ibin=1; ibin<=NBINS_E; ibin++){
		double E_up = Eup.GetBinContent(ibin);
		double E_dn = Edn.GetBinContent(ibin);
		
		int idx1 = (ibin-1)*NBINS_ELECTRONIC;
		
		for(int jbin=1; jbin<=NBINS_ELECTRONIC; jbin++){
			double t = electronic_up.GetBinCenter(jbin);

			double weight_up = E_up*pulse_shape_pre[jbin+idx1];
			double weight_dn = E_dn*pulse_shape_pre[jbin+idx1];
			electronic_up.Fill(t, weight_up);
			electronic_dn.Fill(t, weight_dn);
		}
	}

	// Calculate fADC values
	double Qup = GetQ(electronic_up);
	double Qdn = GetQ(electronic_dn);
	
	fADC_up = (int)(Qup/0.100); // 100fC LSB for CAEN module (n.b for 4096 counts/2V, this will be 2.5 times larger!)
	fADC_dn = (int)(Qdn/0.100); // 100fC LSB for CAEN module (n.b for 4096 counts/2V, this will be 2.5 times larger!)

	geometric_mean = sqrt((double)fADC_up * (double)fADC_dn);
}

//---------------
// GetQ
//---------------
double DBCALTimeSpectrum::GetQ(DHistogram &h_electronic)
{
	// Calculate signal size in pC from pulse shape
	
	// integral in mV * ns. Note we must multiply by the bin width explicitly!
	double integral = h_electronic.Integral()*h_electronic.GetBinWidth();

	// Convert to nC by dividing by 1000 (for the mV)
	// and dividing by 50 Ohms. (don't divide away
	// the nano in ns so that it stays in nC)
	double Q = integral/1000.0/50.0; // 50 Ohms

	// Signal shape of histogram is for TDC. We want this to
	// represent charge seen by fADC so divide by 10 to 
	// account for different gains.
	Q /=10.0;
	
	// Scale from nC to pC
	Q *= 1000.0;
	
	return Q;
}

//---------------
// MakeTSpline
//---------------
TSpline* DBCALTimeSpectrum::MakeTSpline(void)
{
	// These points are taken by mapping out figure 6
	// of GlueX-doc-1795 (the "5ns" rise time pulse).

	double t[] ={
		-7.00,
		-4.92,
		-0.75,
		4.24,
		9.37,		// 5
		
		9.71,
		10.05,
		10.39,
		10.73,
		
		11.07,
		12.03,
		12.92,
		14.08,
		16.05,	// 10
		18.56,
		21.22,
		24.31,
		28.44,
		32.92,	// 15
		38.64,
		42.07,
		//47.10,
		//52.59,
		//58.22,	// 20
		67.95,
		71.74,
		79.85,
		120.0,
		125.0,	// 25
		130.0,
		160.0,
		180.0,
		-1000.0 // flag end of data
	};
	
	double mV[]={
		0,
		0,
		0,
		0,
		0-1,					// 5

		-1.2,
		-1.4,
		-1.6,
		-1.8,

		0-2,
		-49.16039931,
		-159.0483507,
		-375.9324653,
		-711.3798959,	// 10
		-902.2379167,
		-988.9915625,
		-1020.801233,
		-960.0736806,
		-850.1857292,	// 15
		-694.0291667,
		-601.4919445,
		//-539.29,
		//-419.88,
		//-300.46,			// 20
		-142.53,
		-100.15,
		-61.63,
		-8.0,
		-4.0,				// 25
		-2.0,
		-1.0,
		0.0,
		-1000.0 // flag end of data
	};

	// Find number points
	int Npoints=0;
	while(true){
		if(t[Npoints]==-1000.0 && mV[Npoints]==-1000.0)break;
		Npoints++;
	}
	
	TSpline *s = new TSpline3("spline", t, mV, Npoints);
	
	return s;
}

//---------------
// FindTimes
//---------------
void DBCALTimeSpectrum::FindTimes(double thresh_mV)
{
	int bin_up = electronic_up.FindFirstBinAbove(thresh_mV);
	int bin_dn = electronic_dn.FindFirstBinAbove(thresh_mV);
	
	tup = bin_up>0 ? electronic_up.GetBinCenter(bin_up):-1000.0;
	tdn = bin_dn>0 ? electronic_dn.GetBinCenter(bin_dn):-1000.0;
	
	tup_corrected = tup - (bin_up>0 ? BCAL_TimeWalkUp(fADC_up, layer):0.0);
	tdn_corrected = tdn - (bin_dn>0 ? BCAL_TimeWalkDn(fADC_dn, layer):0.0);
}

//---------------
// Reset
//---------------
void DBCALTimeSpectrum::Reset(void)
{				
	Eup.Reset();
	Edn.Reset();
	electronic_up.Reset();
	electronic_dn.Reset();

	Etot = 0.0;
	fADC_up = 0;
	fADC_dn = 0;
	geometric_mean = 0.0;
	tup = tdn = -1000.0;
	tup_corrected = tdn_corrected = -1000.0;
}

//---------------
// SamplingSmear
//---------------
void DBCALTimeSpectrum::SamplingSmear(double Emax, double theta_thrown)
{
	// Here we apply a systematic error due to the sampling
	// fluctuations. The total energy (Emax) is taken from the
	// highest energy particle in the DMCTrajectoryPoints, but is
	// not used. We calculate a sigma based on the deposited energy
	// and given theta.
	//
	// The error is applied by finding the ratio of the smeared
	// cell energy to unsmeared cell energy and scaling both Eup and Edn
	// by it. Note that this does not affect the electronic 
	// histograms since it is assumed they will be calculated
	// afterwards from Eup and Edn.
	
	// Etot must be converted to GeV
	double Etot_GeV = Etot/1000.0;
	
	// Sampling fluctuations from calibDB
	double BCAL_SAMPLINGCOEFA = 0.042;
	double BCAL_SAMPLINGCOEFB = 0.013;
	double sigmaSamp = BCAL_SAMPLINGCOEFA / sqrt( Etot_GeV ) + BCAL_SAMPLINGCOEFB;

	// From plots in an e-mail sent by Matt S. on July 14th or 15th
	// indicating an angular dependance for 150 MeV photons
	double BCAL_SAMPLING_THETA_A = 0.01046;
	double BCAL_SAMPLING_THETA_B = 0.11260;
	sigmaSamp *= (BCAL_SAMPLING_THETA_A/sin(theta_thrown) + BCAL_SAMPLING_THETA_B)/(BCAL_SAMPLING_THETA_A + BCAL_SAMPLING_THETA_B);

	// The value of sigmaSamp is now actually sigma/E. Multiply by
	// Etot to get it in GeV.
	sigmaSamp *= Etot_GeV;
	
	// Randomly sample the fluctuation
	double Esmeared = rnd.Gaus(Etot_GeV,sigmaSamp);

	// Calculate ratio of smeared to unsmeared
	double ratio = Esmeared/Etot_GeV;

	// Scale attenuated energy histos for each end
	Eup.Scale(ratio);
	Edn.Scale(ratio);
}

//---------------
// PoissonSmear
//---------------
void DBCALTimeSpectrum::PoissonSmear(void)
{
	// Apply Poisson statistics due to photo-electrons
	// to the attenuated energy spectra. We do this by
	// converting the integral of the attenuated energy
	// into photo electrons and then sampling from a
	// Poisson distribution with that mean. The ratio
	// of the quantized, sampled value to the unquantized
	// integral (in PE) is used to scale the spectrum.

	double MeV_per_PE = 0.668; // see slides shown at 7/21/2011 BCAL segmentation meeting

	// Integrate Eup and Edn
	double Iup = Eup.Integral(); // in attenuated MeV
	double Idn = Edn.Integral(); // in attenuated MeV

	if(Iup>0.0){
		// Convert to number of PE
		double mean_pe_up = Iup/MeV_per_PE;
		
		int Npe_up = rnd.Poisson(mean_pe_up);
		double ratio_up = (double)Npe_up/mean_pe_up;
		
		Eup.Scale(ratio_up);
	}

	if(Idn>0.0){
		// Convert to number of PE
		double mean_pe_dn = Idn/MeV_per_PE;
		
		int Npe_dn = rnd.Poisson(mean_pe_dn);
		double ratio_dn = (double)Npe_dn/mean_pe_dn;
		
		Edn.Scale(ratio_dn);
	}
}

//---------------
// DetectorJitterSmear
//---------------
void DBCALTimeSpectrum::DetectorJitterSmear(void)
{
	// The photo detection device has an intrinsic timing resolution
	// caused by jitter in the time between the input signal (light)
	// and the generated electronic pulse. This is quoted in the 
	// spec sheet as "Time Resolution (FWHM)" for a single photon.
	// (see page 2 of the following link:
	// http://sales.hamamatsu.com/assets/pdf/parts_S/s10362-33series_kapd1023e05.pdf)
	//
	// n.b. XP2020 PMT's have a very similar time jitter as the SiPMs
	//
	// To include this effect, each bin is converted into a number
	// of photo-electrons and each of them is randomly displaced in time
	// with the specified width. Because the photo-statistics is
	// applied on a global scale, we do not re-quantize the bin contents
	// into an integer number of photo-electrons. Rather, we split it
	// into a number of pieces that is equal to the approximate number
	// of photo-electrons coming from this bin and shift those pieces.
	// Bins with trace amounts of energy are still kept and shifted
	// as well to keep the energy integral constant.
	
	double fwhm_jitter = 0.600; // ns
	double sigma_jitter = fwhm_jitter/2.35; // ns
	double MeV_per_PE = 0.668; // see slides shown at 7/21/2011 BCAL segmentation meeting

	// Upstream
	DHistogram tmp(NBINS_E, E_LO, E_HI);
	for(int ibin=1; ibin<=NBINS_E; ibin++){
		double E = Eup.GetBinContent(ibin);
		if(E==0.0)continue;
		

		int Npe = 1 + (int)floor(E/MeV_per_PE);
		double dE = E/(double)Npe;
		double tbin = Eup.GetBinCenter(ibin);
		
		for(int i=0; i<Npe; i++){
			double t = rnd.Gaus(tbin, sigma_jitter);
			tmp.Fill(t, dE);
		}
	}
	Eup.Reset();
	Eup.Add(&tmp);

	// Downstream
	tmp.Reset();
	for(int ibin=1; ibin<=NBINS_E; ibin++){
		double E = Edn.GetBinContent(ibin);
		if(E==0.0)continue;
		

		int Npe = 1 + (int)floor(E/MeV_per_PE);
		double dE = E/(double)Npe;
		double tbin = Edn.GetBinCenter(ibin);
		
		for(int i=0; i<Npe; i++){
			double t = rnd.Gaus(tbin, sigma_jitter);
			tmp.Fill(t, dE);
		}
	}
	Edn.Reset();
	Edn.Add(&tmp);

}