########################################################################
# Time of flight start detector
########################################################################
processor Lycca/ToFStart/Preproc   UTILS.MhTdcPreprocessor
	input[0:15]     <-   LyccaTargetTofCrate.mhtdc0[0:15]
	input[16:31]    <-   LyccaTargetTofCrate.mhtdc1[16:31]
	# ...
end
processor Lycca/ToFStart/Membrane         LYCCA.CircularMembraneScintillator
	pmt_time[0:31]   <-    Lycca/ToFStart/Preproc.output[0:31]
	x_hit            <-    Frs/S4tracking.xs[5]
	y_hit            <-    Frs/S4tracking.ys[5]
        #....
end
# create the delata_t vs. distance corellation plots for all 32 PMTs
# Membrane.dist[i] is the distance of PMT_i to the particle impact point
# Membrane.T_Tp[i] is the time of PMT_i minus the estimated time of particle impact
for $i in [0:31]
	processor Lycca/ToFStart/Diagnostics/Tdist0$i     UTILS.Pair
		first   <-   Lycca/ToFStart/Membrane.dist[$i]
		second  <-   Lycca/ToFStart/Membrane.T_Tp[$i]
		display first:second 256,0,300:512,-35,55       in   TofStart_Tdist_corrlation
	end
end

#########################################################################
# The following two processors extract the Sc41 time
# Don't forget to set the gates around the peak in the preprocessing
#########################################################################
processor Lycca/ToFfrs/SciPreproc   UTILS.MhTdcPreprocessor
	input[0]     <-   LyccaTargetTofCrate.mhtdc0[16] // Sci21
	input[1]     <-   LyccaTargetTofCrate.mhtdc0[18] // Sci21
	input[2]     <-   LyccaTargetTofCrate.mhtdc0[20] // Sci41
	input[3]     <-   LyccaTargetTofCrate.mhtdc0[22] // Sci41
	display diff   1,0.1,4096,200  | diff_gate   in Preproc/diff
	display output 1,0.1,4096,200                in Preproc/output
end
processor Lycca/ToFfrs/Sc41  UTILS.Pair
	first  <- Lycca/ToFfrs/SciPreproc.output[2]
	second <- Lycca/ToFfrs/SciPreproc.output[3]
	display average
end
#########################################################################
# The following two processors do the normal ToFStart Membrane processing
#########################################################################
processor Lycca/ToFStart/Preproc   UTILS.MhTdcPreprocessor
	input[0:15]     <-   LyccaTargetTofCrate.mhtdc0[0:15]
	input[16:31]    <-   LyccaTargetTofCrate.mhtdc1[16:31]
	display diff 1,0.1,4096,200  | diff_gate   in Preproc/diff
	display output 1,0.1,4096,200              in Preproc/output
end
processor Lycca/ToFStart/Membrane         LYCCA.CircularMembraneScintillator
	pmt_time[0:31]   <-    Lycca/ToFStart/Preproc.output[0:31]
	x_hit            <-    Frs/S4tracking.xs[5]
	y_hit            <-    Frs/S4tracking.ys[5]
	display   x:y               256,-150,150:256,-150,150
	display x_hit:x
	display y_hit:y
end
########################################################################
# Do the correlation of  Sc41-PMTi vs. PMTdistance and determine the 
# slope of the line the 2d histograms
########################################################################
for $i in [0:31]
	processor Lycca/ToFStart/Diagnostics/TSc41_diff$i       UTILS.Pair
		first   <- Lycca/ToFStart/Membrane.pmt_time[$i]
		second  <- Lycca/ToFfrs/Sc41.average
	end
	processor Lycca/ToFStart/Diagnostics/TSc41_dist$i     UTILS.Pair
		first   <-   Lycca/ToFStart/Membrane.dist[$i]
		second  <-   Lycca/ToFStart/Diagnostics/TSc41_diff$i.difference
		display first:second  100,0,300:200,-15,15  in   TofStart_TSc41_dist_corrlation
	end
end

	if (time_computed == false || parameter(always_compute_position) != 0)
	{
		// We haven't been informed about the position of the particle impact
		// position or we want to compute the particle position in any case. 
		// It will be determined by the following minimum search algorithm.
		double radius_sqr = parameter(radius)*parameter(radius);
		
		// Tiny random jitter around the beam axis is the starting point.
		x_particle = 0.2*parameter(radius)*(rand()%2000 - 1000)/1000.0;
		y_particle = 0.2*parameter(radius)*(rand()%2000 - 1000)/1000.0;
		
		// Select the start step size more or less arbitrary.
		// 1/7 of the disk radius seems to be a reasonable choice.
		double stepsize = parameter(radius)/7;
		int direction = 0;
		double last_variance = analyze(x_particle, y_particle, true);
		for (int n = 0;; ++n)
		{
			if (x_particle*x_particle + y_particle*y_particle > 4*radius_sqr) // we left the circular membrane
			{
				//std::cerr << n << " steps to break the circle at position "
				//		  << x_particle << " " << y_particle
				//		  << "   stepsize was " << stepsize << std::endl;
				return;
			}
			double variance;
			// From the current position go strictly downhill
			// in one of the four directions (+x,-x,+y,-y)
			// on the chi^2 plane.
			if (     (variance = analyze(x_particle + stepsize, y_particle,            true)) < last_variance)
				x_particle += stepsize;
			else if ((variance = analyze(x_particle - stepsize, y_particle,            true)) < last_variance)
				x_particle -= stepsize;
			else if ((variance = analyze(x_particle,            y_particle + stepsize, true)) < last_variance)
				y_particle += stepsize;
			else if ((variance = analyze(x_particle,            y_particle - stepsize, true)) < last_variance)
				y_particle -= stepsize;
			else
			{
				// If a minimum is found and the stepsize is not small
				// enough, reduce the stepsize and restart searching.
				if (stepsize > 2.0)
				{
					stepsize /= 2;
				}
				else
				{
					// If the stepsize was small enough (less than 2mm
					// seems to be reasonable, because the spacial
					// resolution of the device is approx. 6-7mm)
					// we are done.
					//std::cerr << n << " steps to find the position " << x_particle << " " << y_particle << std::endl;
					break;
				}
			}
			
			last_variance = variance;
		}