// $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
static double BCAL_TimeWalkUp(double fADC, int layer);
static double BCAL_TimeWalkDn(double fADC, int layer);
#include "timewalk_FINE.h"

// This is defined in DBCAL_TimeSpectrum_factory.cc.
// It can be changed via configuration parameter
// MEV_PER_PE.
extern double MeV_per_PE;

// This is defined in DBCAL_TimeSpectrum_factory.cc.
// It can be changed via configuration parameter
// APPLY_ELECTRONIC_NOISE.
extern bool APPLY_ELECTRONIC_NOISE;

// This is defined in DBCAL_TimeSpectrum_factory.cc.
// It can be changed via configuration parameter
// ADDITIONAL_GAIN.
extern double ADDITIONAL_GAIN;

// Sampling fluctuation parameters. Can be set via
// configuration parameter.
extern double BCAL_SAMPLINGCOEFA;
extern double BCAL_SAMPLINGCOEFB;

// The preamp for the final BCAL has 2 stages. The first
// has a gain of 2 according to Fernando (conversation
// on 10/10/2011).
double preamp_gain_stage1 = 2.0;

// Set in DBCALTimeSpectrum_factory.cc via config parameter.
extern double PREAMP_GAIN_TDC;

// Conversion factor from charge in pC to fADC counts
// This is based on a single sample(4ns) of the 12-bit 
// fADC250 being 2V maximum:
// I = (2V/50Ohms)  ;  Q = I*4ns  ;  LSB = Q/4096
double LSB_pC_fADC = 2.0/50.0*4.0E3/4096.0;

//---------------
// 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;
	}

	// Optionally save the before and after electronic pulse shapes
	bool save_elec_pulse = false;
	if(save_elec_pulse){
		char hname[256];
		sprintf(hname, "elec_up_before_lay%0d_sec%0d", layer,sector);
		electronic_up.MakeTH1D(hname, "Upstream Before");
		sprintf(hname, "elec_dn_before_lay%0d_sec%0d", layer,sector);
		electronic_dn.MakeTH1D(hname, "Downstream Before");
	}

	// Add in electronic noise due to summing circuit plus cable pickup
	if(APPLY_ELECTRONIC_NOISE){
		double noise_rate_GHz = 1.0;
		double noise_amp_mV = 1.0;
		AddElectronicNoise(noise_rate_GHz, noise_amp_mV, 2.0);
	}

	// Optionally save the before and after electronic pulse shapes
	if(save_elec_pulse){
		char hname[256];
		sprintf(hname, "elec_up_after_lay%0d_sec%0d", layer,sector);
		electronic_up.MakeTH1D(hname, "Upstream After");
		sprintf(hname, "elec_dn_after_lay%0d_sec%0d", layer,sector);
		electronic_dn.MakeTH1D(hname, "Downstream After");
	}

	// 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);
	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;

	// Get primary dark pulses
	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 ));

	// Randomly place photo-electrons in time
	for(int i=0; i<darkPulse_up; i++){
		double s  = rnd.Rndm();
		double t  = tmin + s*BCAL_INTWINDOW_NS;
		double dE = MeV_per_PE;  // single photo electron
		if(rnd.Rndm() <= BCAL_XTALK_FRACT){
			dE += MeV_per_PE; // single cross-talk pixel
			if(rnd.Rndm() <= BCAL_XTALK_FRACT){
				dE += MeV_per_PE; // second cross-talk pixel
			}
		}
		Eup.Fill(t, dE); 
	}
	for(int i=0; i<darkPulse_dn; i++){
		double s  = rnd.Rndm();
		double t  = tmin + s*BCAL_INTWINDOW_NS;
		double dE = MeV_per_PE;  // single photo electron
		if(rnd.Rndm() <= BCAL_XTALK_FRACT){
			dE += MeV_per_PE; // single cross-talk pixel
			if(rnd.Rndm() <= BCAL_XTALK_FRACT){
				dE += MeV_per_PE; // second cross-talk pixel
			}
		}
		Edn.Fill(t, dE); 
	}	
}

//---------------
// 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
			// (see slides from BCAL segmentation meeting 7/21/2011)
			double QCD_counts_per_PE = 61.67;
			double LSB_pC_CAEN_V792 = 0.100;
			double gain_factor_for_test = 67.0;
			double TDC_gain_factor = preamp_gain_stage1*PREAMP_GAIN_TDC*ADDITIONAL_GAIN;
			double SiPM_pulse_integral_pC = 1200.0;
			double SiPM_pulse_peak_mV = 2293.0;
			
			double mV_per_MeV = 1.0; // just a factor to start with so we can multiply/divide easily below
			mV_per_MeV *= QCD_counts_per_PE;			// QCD counts/PE
			mV_per_MeV *= LSB_pC_CAEN_V792;			// pC/PE
			mV_per_MeV /= MeV_per_PE;					// pC/MeV
			mV_per_MeV /= gain_factor_for_test;		// divide out gain of pre-amp Carl/Yi's setup used
			mV_per_MeV *= TDC_gain_factor; 			// multiply total gain of pre-amp for TDC leg
			mV_per_MeV *= SiPM_pulse_peak_mV/SiPM_pulse_integral_pC; // mV/MeV
	
			// Normalize pulse shape. 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 = -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/LSB_pC_fADC);
	fADC_dn = (int)(Qdn/LSB_pC_fADC);

	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 PREAMP_GAIN_TDC to 
	// account for different gains.
	Q /=PREAMP_GAIN_TDC;
	
	// 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)
{
	// Electronic noise will smear out the threshold-signal
	// noise difference. Implement that effect here by
	// just smearing the threshold.
	thresh_mV += rnd.Gaus(0.0, 2.0); // smear by 2mV electronic noise

	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
	// (coming from class member Etot) and theta_thrown.
	//
	// 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.
	//
	// This has been reworked a bit (10/20/2011) to add the floor term
	// and stochastic term in quadrature rather than linearly.
	
	// If no energy is deposited, return now
	if(Etot==0.0)return;
	
	// Etot must be converted to GeV
	double Etot_GeV = Etot/1000.0;

	// Theta dependence of sampling fluctuations
	// 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;
	double sigma_scale = (BCAL_SAMPLING_THETA_A/sin(theta_thrown) + BCAL_SAMPLING_THETA_B)/(BCAL_SAMPLING_THETA_A + BCAL_SAMPLING_THETA_B);

	double sigma_stochastic = BCAL_SAMPLINGCOEFA*sqrt( Etot_GeV );
	sigma_stochastic *= sigma_scale; // include theta dependence
	double dE_stochastic = rnd.Gaus(0.0, sigma_stochastic);

	double sigma_floor = BCAL_SAMPLINGCOEFB*Etot_GeV;
	sigma_floor *= sigma_scale; // include theta dependence
	double dE_floor = rnd.Gaus(0.0, sigma_floor);
	
	double dE = dE_stochastic + dE_floor;
	
	double Esmeared = Etot_GeV + dE;

	// 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.

	// 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

	// 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);

}


//---------------
// AddElectronicNoise
//---------------
void DBCALTimeSpectrum::AddElectronicNoise(double rate_GHz, double amp_sigma_mV, double max_ramp_rate_mV_per_ns)
{
	// This is used to add noise to the electronic_up and electronic_dn
	// histograms that arises from the electronics. Specifically, from
	// the summing circuit as that will affect different segmentation
	// schemes.
	//
	// This works by first defining a spline using points whose distribution
	// in time is exponential. The time constant is given by 1/rate_GHz.
	// Time differences are randomly selected and then an amplitude is
	// randomly selected from a Gaussian whose sigma is given by
	// amp_sigma_mV. The ramp rate from the last point to the new point
	// is calculated (inmV/ns) and if it exceeds max_ramp_rate_mV_per_ns
	// then the point is discarded and a new one chosen. Points with zero
	// amplitude are always added as the first and last points in the
	// spline to make sure it is well-behaved throughout the entire time
	// window.
	// 
	// The spline is then evaluated for each bin center in the histogram
	// to get the noise for that bin which is added to the bin. Separate
	// splines are created for the up and dn histos.

	double t0 = ELECTRONIC_LO;
	double t_width = ELECTRONIC_HI - ELECTRONIC_LO;
	
	TSpline *sp_up = MakeEnoiseTSpline(t_width, rate_GHz, amp_sigma_mV, max_ramp_rate_mV_per_ns);
	TSpline *sp_dn = MakeEnoiseTSpline(t_width, rate_GHz, amp_sigma_mV, max_ramp_rate_mV_per_ns);

	// Loop over all bins of electronic histos, adding electronic noise
	for(int ibin=1; ibin<=NBINS_ELECTRONIC; ibin++){
		double c_up = electronic_up.GetBinContent(ibin);
		double c_dn = electronic_dn.GetBinContent(ibin);

		double t = electronic_up.GetBinCenter(ibin) - t0;

		c_up += sp_up->Eval(t);
		c_dn += sp_dn->Eval(t);
		
		electronic_up.SetBinContent(ibin, c_up);
		electronic_dn.SetBinContent(ibin, c_dn);
	}

	delete sp_up;
	delete sp_dn;
}

//---------------
// MakeEnoiseTSpline
//---------------
TSpline* DBCALTimeSpectrum::MakeEnoiseTSpline(double t_width, double rate_GHz, double amp_sigma_mV, double max_ramp_rate_mV_per_ns)
{
	// This creates a TSpline representing electronic noise 
	// defined by the input parameters. See comments at top of
	// AddElectronicNoise() for more.

	// Start off at t=0 and randomly choose a time difference
	// to the next point assuming an exponential time structure
	// until we get past the end of the time window.
	double t = 0.0;
	vector<double> times;
	vector<double> amps;
	times.push_back(0); // add first point as 0 amplitude
	amps.push_back(0.0);
	do{
		double s = rnd.Rndm();
		double delta_t = -1.0/rate_GHz * log(1.0-s);
		t += delta_t;
		if(t<t_width){
			double last_t = times[amps.size()-1];
			double last_amp = amps[amps.size()-1];

			double A = rnd.Gaus(0.0, amp_sigma_mV);
			
			// Calculate rise time and check that it is not too steep
			double ramp_rate = (A-last_amp)/(t - last_t);
			if(fabs(ramp_rate) <= max_ramp_rate_mV_per_ns){
				times.push_back(t);
				amps.push_back(A);
			}else{
				// Ramp rate is too high. Discard this time
				// and try again next iteration
				t -= delta_t;
			}
		}
		
	}while(t<t_width);
	times.push_back(t_width); // add last point as 0 amplitude
	amps.push_back(0.0);
	
	TSpline3 *spline = new TSpline3("E-noise", &times[0], &amps[0], times.size());

	return spline;
}