#include "DSSSD.hpp"
#include <cmath>
#include <iostream>

#include "fit.hpp"


DSSSD::DSSSD(const std::string &config_dir,
		const std::string &name,
		bool no_files)
	: Processor(name, config_dir, no_files)
{
	// setup of parameters we expect to find in the parameter file for DSSSD
	NAME_PARAMETER(strip_num_n);
	NAME_PARAMETER(thresholds_n);
	NAME_PARAMETER(strip_num_p);
	NAME_PARAMETER(thresholds_p);
	NAME_PARAMETER(x_center);
	NAME_PARAMETER(y_center);
	NAME_PARAMETER(strip_width_n);
	NAME_PARAMETER(strip_width_p);
	NAME_PARAMETER(left_right);
	NAME_PARAMETER(top_bottom);
	NAME_PARAMETER(rotate);
	NAME_PARAMETER(selfcalibration);
	NAME_PARAMETER(selfcalibration2);
	// read the parameter file before setting the input/output channels
	// because they depend on some of the parameters
	init();
	read_parameters();

	// specify in and output
	NAME_INPUT_ARRAY(amplitude_p, parameter(strip_num_n));
	NAME_INPUT_ARRAY(amplitude_n, parameter(strip_num_p));
	NAME_INPUT_ARRAY(time_p,      parameter(strip_num_p));

	NAME_OUTPUT_ARRAY(cal_amplitude_p, parameter(strip_num_n));
	NAME_OUTPUT_ARRAY(cal_amplitude_n, parameter(strip_num_p));
	NAME_OUTPUT_ARRAY(cal_time_p,      parameter(strip_num_p));

	NAME_OUTPUT_CHANNEL(x);
	NAME_OUTPUT_CHANNEL(y);
	NAME_OUTPUT_CHANNEL(x_sub);
	NAME_OUTPUT_CHANNEL(y_sub);
	NAME_OUTPUT_CHANNEL(dE);
	NAME_OUTPUT_CHANNEL(dE_p);
	NAME_OUTPUT_CHANNEL(dE_n);
	NAME_OUTPUT_CHANNEL(multiplicity_p);
	NAME_OUTPUT_CHANNEL(multiplicity_n);
	NAME_OUTPUT_CHANNEL(cluster_multiplicity_p);
	NAME_OUTPUT_CHANNEL(cluster_multiplicity_n);
	NAME_OUTPUT_CHANNEL(dE_max_p);		// Added by LS
	NAME_OUTPUT_CHANNEL(dE_max_n);		// Added by LS
	NAME_OUTPUT_CHANNEL(dE_max);		// Added by LS
// 	NAME_OUTPUT_CHANNEL(dE_max_cal_p);	// Added by LS
// 	NAME_OUTPUT_CHANNEL(dE_max_cal_n);	// Added by LS

	// setup calibration
	NAME_COEFFICIENTS_ARRAY(cal_p, parameter(strip_num_p));
	NAME_COEFFICIENTS_ARRAY(cal_n, parameter(strip_num_n));
	NAME_COEFFICIENTS_ARRAY(cal_t, parameter(strip_num_p));

	NAME_COEFFICIENTS(cal_dE);
	
	init();
	read_parameters();
	read_coefficients();
	
	if (parameter_count(selfcalibration) == 4 && parameter(selfcalibration, 0) > 0)
	{
		const int N = parameter(selfcalibration, 0);
		prob_logs_np.resize(parameter(strip_num_n),
		                    std::vector<std::vector<double> >(parameter(strip_num_p), 
		                                                      std::vector<double>(N, 1.0/N) ) );
		                                                      
		slopes.resize(parameter(strip_num_n), std::vector<double>(parameter(strip_num_p), 1));
		slope_errors.resize(parameter(strip_num_n), std::vector<double>(parameter(strip_num_p), 1));
	}
	if (parameter_count(selfcalibration2) == 4 && parameter(selfcalibration2, 0) > 0)
	{
		const int N = parameter(selfcalibration2, 0);
		prob_log_n_slopes.resize(parameter(strip_num_n), std::vector<double>(N, 1.0/N));
		prob_log_p_slopes.resize(parameter(strip_num_p), std::vector<double>(N, 1.0/N));
	}
}

DSSSD::~DSSSD()
{
}


void DSSSD::process(prespec::viscon::Interface &viscon_interface, int trigger)
{
	if (input_array_size(amplitude_n) == 0 ||
		input_array_size(amplitude_p) == 0)
		return;

	// Do a selfconsitent automatic selfcalibration.
	// This will create a file with coefficients that can be directly used
	// as calibration for this processor
	if (parameter_count(selfcalibration) == 4 && parameter(selfcalibration, 0) > 0)
	{
		do_selfcalibration();
	}

	// treat the time signals
	for (int i = 0; i < input_array_size(time_p); ++i)
	{
		int index = input_array_index(time_p, i);
		double t_raw = input_array_value(time_p, i);
		double t_cal = calibrate(t_raw, cal_t, index);
		fill_output_array(cal_time_p, index, t_cal);
	}

	int y_factor = (parameter(top_bottom)>=0)?1:-1;
	int x_factor = (parameter(left_right)>=0)?1:-1;
	// process the p-side
	// variables for maximum-energy-deposition  position determination
	int num_indices_p = 0;
	int max_index_p = -1;
	double max_dE_p = 0;
	// variables for energy-weighted-sub-pixel  position determination
	double mean_index_p = 0;
	double sum_dE_p = 0;
	int mult_p = 0;
	std::vector<int> hitpat_p(parameter(strip_num_p),0);
	std::vector<double> dEraw_p(parameter(strip_num_p),0);
	for (int i = 0; i < input_array_size(amplitude_p); ++i)
	{
		int index = input_array_index(amplitude_p, i);
		double dEraw = input_array_value(amplitude_p, i);

		if (dEraw <= parameter(thresholds_p,0) ||
			dEraw >= parameter(thresholds_p,1))	continue;

		++mult_p;
		hitpat_p[index] = 1;

		double dE = calibrate(amplitude_p, index, dEraw);
		dEraw_p[index] = dE;
		fill_output_array(cal_amplitude_p, index, dE);
		++num_indices_p;
		mean_index_p += index * dE;
		sum_dE_p += dE;
		if (i == 0 || max_dE_p < dE)
		{
			max_dE_p = dE;
			max_index_p = index;
		}
	}
	set_output(multiplicity_p, mult_p);
	
	// =======  Setting up max E sums for sub-pixel algorithm	=======	// Added by LS
	double mean_max_index_p = 0;
	double sum_max_dE_p = 0;
	int max_indices_p = 0;
	if( max_index_p > 0 && max_index_p < (parameter(strip_num_p)-1) )
	{
	  for (int i=(max_index_p-1); i<(max_index_p+2); i++)
	  {
		mean_max_index_p += i * dEraw_p[i];
		sum_max_dE_p += dEraw_p[i];
		if( dEraw_p[i] > 0 )
		{
			++max_indices_p;
		}
	  }
	}
	if( max_index_p == 0 )
	{
	  for (int i=max_index_p; i<(max_index_p+2); i++)
	  {
		mean_max_index_p += i * dEraw_p[i];
		sum_max_dE_p += dEraw_p[i];
		if( dEraw_p[i] > 0 )
		{
			++max_indices_p;
		}
	  }
	}
	if( max_index_p == (parameter(strip_num_p)-1) )
	{
	  for (int i=(max_index_p-1); i<(max_index_p+1); i++)
	  {
		mean_max_index_p += i * dEraw_p[i];
		sum_max_dE_p += dEraw_p[i];
		if( dEraw_p[i] > 0 )
		{
			++max_indices_p;
		}
	  }
	}
	// =======   End of addition   =======
	
	// process the n-side
	// variables for maximum-energy-deposition  position determination
	int num_indices_n = 0;
	int max_index_n = -1; // for strip resolution
	double max_dE_n = 0;  // (get strip with maximum energy deposition)
	// vairables for energy-weighted-sub-pixel  position determination
	double mean_index_n = 0; // for sub strip resolution
	double sum_dE_n = 0;     // (calculate enery weighted mean position)
	int mult_n = 0;
	std::vector<int> hitpat_n(parameter(strip_num_n),0);
	std::vector<double> dEraw_n(parameter(strip_num_n),0);
	for (int i = 0; i < input_array_size(amplitude_n); ++i)
	{
		int index = input_array_index(amplitude_n, i);
		double dEraw = input_array_value(amplitude_n, i);

		if (dEraw <= parameter(thresholds_n,0) ||
			dEraw >= parameter(thresholds_n,1))	continue;

		++mult_n;
		hitpat_n[index] = 1;
		dEraw_n[index] = dE;

		double dE = calibrate(amplitude_n, index, dEraw);
		fill_output_array(cal_amplitude_n, index, dE);
		dEraw_n[index] = dE;
		++num_indices_n;
		
		mean_index_n += index * dE;
		sum_dE_n += dE;
		if (i == 0 || max_dE_n < dE)
		{
			max_dE_n = dE;
			max_index_n = index;
		}
	}
	set_output(multiplicity_n, mult_n);

	
	//  ======    Setting up max E sums for sub-pixel algorithm  ========	// Added by LS
	double mean_max_index_n = 0;
	double sum_max_dE_n = 0;
	int max_indices_n = 0;
	if( max_index_n > 0 && max_index_n < (parameter(strip_num_n) - 1))
	{
	  for (int i=(max_index_n-1); i<(max_index_n+2); i++)
	  {
		mean_max_index_n += i * dEraw_n[i];
		sum_max_dE_n += dEraw_n[i];
		if( dEraw_n[i] > 0 )
		{
			++max_indices_n;
		}
	  }
	}
	
	if( max_index_n == 0 )
	{
	  for (int i=max_index_n; i<(max_index_n+2); i++)
	  {
		mean_max_index_n += i * dEraw_n[i];
		sum_max_dE_n += dEraw_n[i];
		if( dEraw_n[i] > 0 )
		{
			++max_indices_n;
		}
	  }
	}
	
	if( max_index_n == (parameter(strip_num_n) - 1) )
	{
	  for (int i=(max_index_n-1); i<(max_index_n+1); i++)
	  {
		mean_max_index_n += i * dEraw_n[i];
		sum_max_dE_n += dEraw_n[i];
		if( dEraw_n[i] > 0 )
		{
			++max_indices_n;
		}
	  }
	}
	// =======   End of addition   =========
	
	int n_clusters_p = 0;
	bool on = false;
	for (int i = 0; i < hitpat_p.size(); ++i)
	{
		if (hitpat_p[i] && !on)
		{
			++n_clusters_p;
			on = true;
		}
		if (!hitpat_p[i])
			on = false;
	}
	set_output(cluster_multiplicity_p, n_clusters_p);

	int n_clusters_n = 0;
	on = false;
	for (int i = 0; i < parameter(strip_num_n); ++i)
	{
		if (hitpat_n[i] && !on)
		{
			++n_clusters_n;
			on = true;
		}
		if (!hitpat_n[i])
			on = false;
	}
	set_output(cluster_multiplicity_n, n_clusters_n);

// 	if (n_clusters_p != 1 || n_clusters_n != 1)	// Original
// 		return;

	// the noninterpolated position
	if (max_dE_n >= 0 && max_dE_p >= 0)
	{
		double x_jitter = static_cast<double>(rand())/RAND_MAX;
		double y_jitter = static_cast<double>(rand())/RAND_MAX;
		double x_pos = x_factor * (x_jitter + max_index_n - parameter(strip_num_n)/2.) * parameter(strip_width_n);
		double y_pos = y_factor * (y_jitter + max_index_p - parameter(strip_num_p)/2.) * parameter(strip_width_p);
		rotate_xy(x_pos, y_pos, x_pos, y_pos, parameter(rotate));
		x_pos += parameter(x_center);
		y_pos += parameter(y_center);
		set_output(x, x_pos);
		set_output(y, y_pos);
	}

	// the interpolated position and energy loss
 	if (sum_max_dE_p > 0 && sum_max_dE_n > 0)  // Added by LS
// 	if (sum_dE_p > 0 && sum_dE_n > 0)			// Original Code

	{
		//Addition by LS for max dE version of measuring x and y
		double x_jitter = (max_indices_n>1)?(0):(static_cast<double>(rand())/RAND_MAX);
		double y_jitter = (max_indices_p>1)?(0):(static_cast<double>(rand())/RAND_MAX);
		double x_pos_sub = x_factor * (x_jitter + mean_max_index_n/sum_max_dE_n - parameter(strip_num_n)/2.) * parameter(strip_width_n);
		double y_pos_sub = y_factor * (y_jitter + mean_max_index_p/sum_max_dE_p - parameter(strip_num_p)/2.) * parameter(strip_width_p);
		// End of addition
		
		// Original Code:
// 		double x_jitter = (num_indices_n>1)?(0):(static_cast<double>(rand())/RAND_MAX);
// 		double y_jitter = (num_indices_p>1)?(0):(static_cast<double>(rand())/RAND_MAX);
// 		double x_pos_sub = x_factor * (x_jitter + mean_index_n/sum_dE_n - parameter(strip_num_n)/2.) * parameter(strip_width_n);
// 		double y_pos_sub = y_factor * (y_jitter + mean_index_p/sum_dE_p - parameter(strip_num_p)/2.) * parameter(strip_width_p);
		// End of Original code
		
		rotate_xy(x_pos_sub, y_pos_sub, x_pos_sub, y_pos_sub, parameter(rotate));
		x_pos_sub += parameter(x_center);
		y_pos_sub += parameter(y_center);
		set_output(x_sub, x_pos_sub);
		set_output(y_sub, y_pos_sub);
		
		// Addition by LS for max dE measurements
		set_output(dE_max_n, sum_max_dE_n);
		set_output(dE_max_p, sum_max_dE_p);
		set_output(dE_max, calibrate(0.5*(sum_max_dE_n+sum_max_dE_p), cal_dE));
		
		// Original Code
		set_output(dE_n, sum_dE_n);
		set_output(dE_p, sum_dE_p);
		set_output(dE, calibrate(0.5*(sum_dE_n+sum_dE_p), cal_dE));
	}
	
// 	if( max_index_n == 0 && max_index_p > 0 )
// 	{
// 		x_pos_sub = x_factor * (x_jitter + max_index_n - parameter(strip_num_n)/2.) * parameter(strip_width_n);
// 		x_pos_sub += parameter(x_center);
// 		set_output(x_sub, x_pos_sub);
// 	}
// 	
// 	if( max_index_p == 0 && max_index_n > 0 )
// 	{
// 		y_pos_sub = y_factor * (y_jitter + max_index_p - parameter(strip_num_p)/2.) * parameter(strip_width_p);
// 		y_pos_sub += parameter(y_center);
// 		set_output(y_sub, y_pos_sub);
// 	}

	// selfcalibration2 requires the dE_p and dE_n outputs
	// so it is done after the DSSSD processing
	if (parameter_count(selfcalibration2) == 4 && parameter(selfcalibration2, 0) > 0)
	{
		do_selfcalibration2();
	}

}

void DSSSD::rotate_xy(double x, double y, double &x_rot, double &y_rot, int i) // i = 0,1,2,3 rotates by phi=0,90,180,270 degrees
{
	switch(i)
	{
		case 0:
			x_rot = x;
			y_rot = y;
		return;
		case 1:
			x_rot = y;
			y_rot = -x;
		return;
		case 2:
			x_rot = -x;
			y_rot = -y;
		return;
		case 3:
			x_rot = -y;
			y_rot = x;
		return;
	}
}

void DSSSD::do_selfcalibration()
{
		static int NUM = 0;
		++NUM;

		// Compute the probability distribution for the np-slope matrix
		// The np-slope matrix has in the n-th line and p-th column
		// the slope of the correlation of the amplitude of the n-th 
		// n-side strip vs the p-th p-side strip.	
		// The probability distribution is computed, using Bayes' theorem.	
		const int    N     = parameter(selfcalibration, 0);
		const double s1    = parameter(selfcalibration, 1);
		const double s2    = parameter(selfcalibration, 2);
		const double width = parameter(selfcalibration, 3);

		for (int n = 0; n < input_array_size(amplitude_n); ++n)
		{
			for (int p = 0; p < input_array_size(amplitude_p); ++p)
			{
				int    n_index = input_array_index(amplitude_n, n);
				double n_value = input_array_value(amplitude_n, n);
				int    p_index = input_array_index(amplitude_p, p);
				double p_value = input_array_value(amplitude_p, p);
				
				if (p_value <= 0 || n_value <= 0)
					continue;
				
				double slope = p_value/n_value;
				double log_slope = log(slope);
				//std::cout << "slope = " << slope << "    log_slope = " << log_slope << std::endl;
				
				double sum = 0;
				for (int i = 0; i < N; ++i)
				{
					double logs = log(s1) + i*(log(s2)-log(s1))/N;
					prob_logs_np[n_index][p_index][i] *= cauchy(logs - log_slope, width);
					sum += prob_logs_np[n_index][p_index][i];
				}
				for (int i = 0; i < N; ++i)
				{
					prob_logs_np[n_index][p_index][i] /= sum;
				}
			}
		}
		
		// Regularly, determine the calibration factors from the pn-slope 
		// matrix and write the resulting coefficients into a file
		if (NUM%10000 == 0)
		{
			// Write the posterior of all the pn-slopes into a file.
			// This is a huge file and it is only useful for demonstation of 
			// how Bayes' theorem works. 
			//std::string name2 = name();
			//for (int i = 0; i < name2.size(); ++i) if (name2[i] == '/') name2[i] = '_';
			//name2.append(".selfcal_probs");
			//std::ofstream out(name2.c_str());
			//for (int i = 0; i < N; ++i)
			//{
			//	double logs = log(s1) + i*(log(s2)-log(s1))/N;
			//	out << exp(logs) << " ";
			//	for (int p = 0; p < parameter(strip_num_p); ++p)
			//	for (int n = 0; n < parameter(strip_num_n); ++n)
			//	{
			//		out << prob_logs_np[n][p][i] << " ";
			//	}
			//	out << std::endl;
			//}
			
			
			// determine the mean value and the variance of the probability distribution
			// of the slope parameter for all the combinations of p- and n-side strips.
			for (int n = 0; n < parameter(strip_num_n); ++n)
			{
				for (int p = 0; p < parameter(strip_num_p); ++p)
				{
					double max_logs = 0;
					double max = 0;	
					double moment0 = 0;
					double moment1 = 0;
					double moment2 = 0;	
					for (int i = 0; i < N; ++i)
					{
						double logs = log(s1) + i*(log(s2)-log(s1))/N;
						double logs2 = log(s1) + (i+1)*(log(s2)-log(s1))/N;
						double s = exp(logs);
						double deltas = exp(logs2) - exp(logs);
						moment0 +=       prob_logs_np[n][p][i] * deltas;
						moment1 += s *   prob_logs_np[n][p][i] * deltas;
						moment2 += s*s * prob_logs_np[n][p][i] * deltas;
						if (i == 0 || max < prob_logs_np[n][p][i])
						{
							max_logs = logs;
							max = prob_logs_np[n][p][i];
						}
					}
					double mean = moment1/moment0;
					double var = moment2/moment0 - mean*mean;
					//out << n << " " << p << " " 
					//    << std::setw(15) << mean << " " 
					//    << std::setw(15) << sqrt(var) << " " 
					//    << std::setw(15) << exp(max_logs) 
					//    << std::endl;
					    
					slopes[n][p]       = mean; 
					slope_errors[n][p] = sqrt(var);   
				}
				//out << std::endl;
			}
			//out << std::endl;
			
			fit_slopes();
		}
}


void DSSSD::fit_slopes()
{
	using namespace gsl::multifit_nlin;
	
	std::vector<Dp<int> > data;
	for (int n = 0; n < parameter(strip_num_n); ++n)
	{
		for (int p = 0; p < parameter(strip_num_p); ++p)
		{
			int coordinate = n*parameter(strip_num_p)+p;
			double value = slopes[n][p];
			double sigma = slope_errors[n][p];
			data.push_back(Dp<int>(coordinate,value,sigma));
		}
	}

	// fill start parameters in array
	// the parameters are the calibration factors (slopes)
	// (the offsets are assumed to be all equal to 0)
	std::vector<double> params;
	for (int i = 0; i < parameter(strip_num_n)+parameter(strip_num_p)-1; ++i)
		params.push_back(1.0);

	// define function to fit and the type of residues to use
	FitFunction function(parameter(strip_num_n), parameter(strip_num_p));
	FitResidues residues;
	bool verbose = false;

	// create a Fit obejct
	Fit<FitFunction, FitResidues, int, std::vector<Dp<int> >::iterator > 
		fit(function, residues, data.begin(), data.end(), &(*params.begin()), &(*params.end()) , verbose);

	double precision = 1e-5;
	int max_number_of_steps = 50;

	// run the fit
	fit.run(precision, max_number_of_steps);
	
	// read the resulting parameters
	fit.get_result_params(params);

	// write the results into a file
	std::string name2 = name();
	for (int i = 0; i < name2.size(); ++i) if (name2[i] == '/') name2[i] = '_';
	name2.append(".selfcalibrate");
	std::ofstream out(name2.c_str());
	for (int i = 0; i < parameter(strip_num_n); ++i)
		out << "cal_n[" << i << "]\t\t0.0 \t" << params[i] << std::endl;
	out << std::endl;	
	for (int i = 0; i < parameter(strip_num_p)-1; ++i)
		out << "cal_p[" << i << "]\t\t0.0 \t" << params[parameter(strip_num_n)+i] << std::endl;
	out << "cal_p[" << parameter(strip_num_p)-1 << "]\t\t0.0 \t1.0" << std::endl;
}


double DSSSD::cauchy(double x, double width)
{
	return width / ( (width*width + x*x) * M_PI );
}

double FitFunction::operator()(int np, const double *pbegin, const double *pend)
{
	// parameter array [pbegin,pend[ contains:
	// sn[0], sn[1], ... , sn[num_n-1], sp[0], sp[1], ... , sp[num_p-2]
	// one parameter (the last sp[num_p-1]) is fixed to 1
	
	int n = np / num_p_;
	int p = np % num_p_;

	double sn = pbegin[n];
	double sp = 1.0; 
	if (n < num_n_-1 && 
	    p < num_p_-1)
	{    
		sp = pbegin[num_n_+p];
	}
	
	// the relation between the slopes and the pn-slope matrix is
	// Snp = sn/sp
	return sn/sp;
}

// use a normal chisquare fit
double FitResidues::operator()(double f, double v, double s)
{
	return (f - v) / s;
}


void DSSSD::do_selfcalibration2()
{
		static int NUM = 0;
		++NUM;

		// Compute the probability distribution for the n-slopes and p-slopes,
		// assuming that the do_selcalibration procedure was already done
		// The probability distribution is computed, using Bayes' theorem.	
		const int    N     = parameter(selfcalibration2, 0);
		const double s1    = parameter(selfcalibration2, 1);
		const double s2    = parameter(selfcalibration2, 2);
		const double width = parameter(selfcalibration2, 3);

		if (output_valid(dE_p))
		{
			for (int n = 0; n < input_array_size(amplitude_n); ++n)
			{
				int    n_index = input_array_index(amplitude_n, n);
				double n_value = input_array_value(amplitude_n, n);
				double energy_p = output_value(dE_p);
				double slope = energy_p / n_value;
				double log_slope = log(slope);

				// update the likelihood function
				double sum = 0;
				for (int i = 0; i < N; ++i)
				{
					double logs = log(s1) + i*(log(s2)-log(s1))/N;
					prob_log_n_slopes[n_index][i] *= cauchy(logs - log_slope, width);
					sum += prob_log_n_slopes[n_index][i];
				}
				// normalize the likelihood function
				for (int i = 0; i < N; ++i)
				{
					prob_log_n_slopes[n_index][i] /= sum;
				}				
			}
		}
		if (output_valid(dE_n))
		{
			for (int p = 0; p < input_array_size(amplitude_p); ++p)
			{
				int    p_index = input_array_index(amplitude_p, p);
				double p_value = input_array_value(amplitude_p, p);
				double energy_n = output_value(dE_n);
				double slope = energy_n / p_value;
				double log_slope = log(slope);

				// update the likelihood function
				double sum = 0;
				for (int i = 0; i < N; ++i)
				{
					double logs = log(s1) + i*(log(s2)-log(s1))/N;
					prob_log_p_slopes[p_index][i] *= cauchy(logs - log_slope, width);
					sum += prob_log_p_slopes[p_index][i];
				}
				// normalize the likelihood function
				for (int i = 0; i < N; ++i)
				{
					prob_log_p_slopes[p_index][i] /= sum;
				}				
			}
		}
			
		// Regularly, determine the calibration factors from the pn-slope 
		// matrix and write the resulting coefficients into a file
		if (NUM%10000 == 0)
		{
			std::vector<double> n_slopes(parameter(strip_num_n));
			std::vector<double> n_slope_errors(parameter(strip_num_n));
			for (int n = 0; n < parameter(strip_num_n); ++n)
			{
				double max_logs = 0;
				double max = 0;	
				double moment0 = 0;
				double moment1 = 0;
				double moment2 = 0;	
				for (int i = 0; i < N; ++i)
				{
					double logs = log(s1) + i*(log(s2)-log(s1))/N;
					double logs2 = log(s1) + (i+1)*(log(s2)-log(s1))/N;
					double s = exp(logs);
					double deltas = exp(logs2) - exp(logs);
					moment0 +=       prob_log_n_slopes[n][i] * deltas;
					moment1 += s *   prob_log_n_slopes[n][i] * deltas;
					moment2 += s*s * prob_log_n_slopes[n][i] * deltas;
					if (i == 0 || max < prob_log_n_slopes[n][i])
					{
						max_logs = logs;
						max = prob_log_n_slopes[n][i];
					}
				}
				double mean = moment1/moment0;
				double var = moment2/moment0 - mean*mean;
					
				n_slopes[n]       = exp(max_logs); 
				n_slope_errors[n] = sqrt(var);   				
			}
			
			std::vector<double> p_slopes(parameter(strip_num_p));
			std::vector<double> p_slope_errors(parameter(strip_num_p));
			for (int p = 0; p < parameter(strip_num_p); ++p)
			{
				double max_logs = 0;
				double max = 0;	
				double moment0 = 0;
				double moment1 = 0;
				double moment2 = 0;	
				for (int i = 0; i < N; ++i)
				{
					double logs = log(s1) + i*(log(s2)-log(s1))/N;
					double logs2 = log(s1) + (i+1)*(log(s2)-log(s1))/N;
					double s = exp(logs);
					double deltas = exp(logs2) - exp(logs);
					moment0 +=       prob_log_p_slopes[p][i] * deltas;
					moment1 += s *   prob_log_p_slopes[p][i] * deltas;
					moment2 += s*s * prob_log_p_slopes[p][i] * deltas;
					if (i == 0 || max < prob_log_p_slopes[p][i])
					{
						max_logs = logs;
						max = prob_log_p_slopes[p][i];
					}
				}
				double mean = moment1/moment0;
				double var = moment2/moment0 - mean*mean;
					
				p_slopes[p]       = exp(max_logs); 
				p_slope_errors[p] = sqrt(var);   				
			}

			// write the results into a file
			std::string name2 = name();
			for (int i = 0; i < name2.size(); ++i) if (name2[i] == '/') name2[i] = '_';
			name2.append(".selfcalibrate2");
			std::ofstream out(name2.c_str());
			for (int i = 0; i < parameter(strip_num_n); ++i)
				out << "cal_n[" << i << "]\t\t0.0 \t" << n_slopes[i] << std::endl;
			out << std::endl;	
			for (int i = 0; i < parameter(strip_num_p); ++i)
				out << "cal_p[" << i << "]\t\t0.0 \t" << p_slopes[i] << std::endl;
			out << std::endl;	
		}	
}