if [1].EQ.'?' .OR. [1].EQ.' ' then
	mess '==============================================================================='
	mess 'FIT_PREC IDLeft B  fromFit toFit t0' 
	mess
	mess 'Kumac to fit the GEANT simulated decay spectra of L,T,R,B.'
	mess 'If IDleft is given, then it is assumed that'
	mess '   IDtop    = IDleft+2,'
	mess '   IDright  = IDleft+3,'
	mess '   IDbottom = IDleft+4.'
	mess
	mess 'B        : magnetic field in gauss'
	mess 'fromFit  : lower boundary for fit'
	mess 'toFit    : upper boundary for fit'
	mess 't0       : time zero'
	mess
	mess '==============================================================================='
	exitm
endif
*
*===========================================================================
current = $hcdir()

fitopt = 'bl0'  | fit option b=parameters with boundaries
                |            l=log likelihood fit (Poisson stat.)
                |            0=do not plot fit result

* define Histo ID's

	idl = [1]
	idt = [idl]+2
	idr = [idl]+3
	idb = [idl]+4
*
*	check if histograms exist in PAW memory
*
	cd //pawc
	
	if ( $hexist([idl]) .eq. 0 ) then
		mess
		mess Did not found Left histogram ID [idl]
		mess
		exitm
	endif 
	if ( $hexist([idt]) .eq. 0 ) then
		mess
		mess Did not found Top histogram ID [idt]
		mess
		exitm
	endif 
	if ( $hexist([idr]) .eq. 0 ) then
		mess
		mess Did not found Right histogram ID [idr]
		mess
		exitm
	endif 
	if ( $hexist([idb]) .eq. 0 ) then
		mess
		mess Did not found Bottom histogram ID  [idb]
		mess
		exitm
	endif 

	shift

* transfer magnetic field strength:

	B = [1]

* get region of histo to be used for fit:

	if ($index([2],'.').NE.0) then
		von = $rsigma([2])
	else
		von = [2]
	endif
	if ($index([3],'.').NE.0) then
		bis = $rsigma([3])
	else
		bis = [3]
	endif

* get rebinning factor and fit options:

	t0 = [4]

	if ( [t0] .eq. ' ' .or. [t0] .eq. '!' ) then
		mess
		mess Missing time zero !
		mess
		exitm
	endif
	if ( $vexist(tzero) ) then
		v/de tzero
	endif
	v/cre tzero(1)
	v/inp tzero(1) [t0]
	        
* check if added histos are to be used and read in histos:

* do the fits:

	freq = [B] * 0.013554	| precession frequency in MHz
	freqmin = 0.8*[freq]
	freqmax = 1.2*[freq]
	freqstep = 0.005*[freq]
*
	if ( $vexist(par) ) then
		v/de par
	endif
	if ( $vexist(pmin) ) then
		v/de pmin
	endif
	if ( $vexist(pmax) ) then
		v/de pmax
	endif
	if ( $vexist(step) ) then
		v/de step
	endif
	if ( $vexist(errpar) ) then
		v/de errpar
	endif
*
*
*
*			  BKG     N0   Asym   mu+ freq.  Phase tau_mu  lambda
*                       =======================================================
*
*	v/cre par(7)  r    0.    500.  0.25       0.        0.  2.19703   0.
*	v/cre pmin(7) r  -10.      0.  0.         0.     -400.  0.        0.
*	v/cre pmax(7) r  100.  10000.  0.60       0.     +400.  3.      100.
*	v/cre step(7) r    0.     10.  0.01       0.        5.  0.        0.  
*	v/cre errpar(7) r  1.      1.  1.         1.        1.  1.        1.
	v/cre par(6)  r    0.    500.  0.25    [freq]      90.  2.19703   
	v/cre pmin(6) r  -10.      0.  0.   [freqmin]    -400.  0.       
	v/cre pmax(6) r  100.  10000.  0.40 [freqmax]    +400.  3.      
	v/cre step(6) r    0.     10.  0.01 [freqstep]      5.  0.         
	v/cre errpar(6) r  0.      1.  1.         1.        1.  0.        

*
	tdec0 = tzero(1)	| vector tzero filled in FIT_PREC.KUMAC
	ttt = Fit from [von] to [bis], t0  = [tdec0]
*	title [ttt]  
*
*
	hi/fit [idl]([von]:[bis]) fit_prec.for [fitopt] 6 par step pmin pmax errpar
*
*	store fit result and error for A_Left, frequency and relaxation rate
*
	a_l    = par(3)
	a_lerr = errpar(3)
	a_l_weight = [a_l]/([a_lerr]*[a_lerr])
	inv_a_lerr = 1./([a_lerr]*[a_lerr])
	f_l    = par(4)
	f_lerr = errpar(4)
	f_l_weight = [f_l]/([f_lerr]*[f_lerr])
	inv_f_lerr = 1./([f_lerr]*[f_lerr])
*	r_l    = par(7)
*	r_lerr = errpar(7)
*	r_l_weight = [r_l]/([r_lerr]*[r_lerr])
*	inv_r_lerr = 1./([r_lerr]*[r_lerr])
*
*	if ( [opt] .eq. 'strexp' ) then
*		beta_l    = par(8)
*		beta_lerr = errpar(8)
*		beta_l_weight = [beta_l]/([beta_lerr]*[beta_lerr])
*		inv_beta_lerr = 1./([beta_lerr]*[beta_lerr])
*	endif
*
*
*	initialize the phase for the top detector
*
	phase = par(5)
	phase = [phase]+90.
	v/inp par(5) [phase]
	hi/fit [idt]([von]:[bis]) fit_prec.for [fitopt] 6 par step pmin pmax errpar
*
*	store fit result and error for A_top
*
	a_t    = par(3)
	a_terr = errpar(3)
	a_t_weight = [a_t]/([a_terr]*[a_terr])
	inv_a_terr = 1./([a_terr]*[a_terr])
	f_t    = par(4)
	f_terr = errpar(4)
	f_t_weight = [f_t]/([f_terr]*[f_terr])
	inv_f_terr = 1./([f_terr]*[f_terr])
*	r_t    = par(7)
*	r_terr = errpar(7)
*	r_t_weight = [r_t]/([r_terr]*[r_terr])
*	inv_r_terr = 1./([r_terr]*[r_terr])
*
*	if ( [opt] .eq. 'strexp' ) then
*		beta_t    = par(8)
*		beta_terr = errpar(8)
*		beta_t_weight = [beta_t]/([beta_terr]*[beta_terr])
*		inv_beta_terr = 1./([beta_terr]*[beta_terr])
*	endif
*
*
*
*	initialize the phase for the right detector
*
	phase = par(5)
	phase = [phase]+90.
	v/inp par(5) [phase]
	hi/fit [idr]([von]:[bis]) fit_prec.for [fitopt] 6 par step pmin pmax errpar

*	store fit result and error for A_right
*
	a_r    = par(3)
	a_rerr = errpar(3)
	a_r_weight = [a_r]/([a_rerr]*[a_rerr])
	inv_a_rerr = 1./([a_rerr]*[a_rerr])
	f_r    = par(4)
	f_rerr = errpar(4)
	f_r_weight = [f_r]/([f_rerr]*[f_rerr])
	inv_f_rerr = 1./([f_rerr]*[f_rerr])
*	r_r    = par(7)
*	r_rerr = errpar(7)
*	r_r_weight = [r_r]/([r_rerr]*[r_rerr])
*	inv_r_rerr = 1./([r_rerr]*[r_rerr])
*
*	if ( [opt] .eq. 'strexp' ) then
*		beta_r    = par(8)
*		beta_rerr = errpar(8)
*		beta_r_weight = [beta_r]/([beta_rerr]*[beta_rerr])
*		inv_beta_rerr = 1./([beta_rerr]*[beta_rerr])
*	endif
*
*
*
*	initialize the phase for the bottom detector
*
	phase = par(5)
	phase = [phase]+90.
	v/inp par(5) [phase]
	hi/fit [idb]([von]:[bis]) fit_prec.for [fitopt] 6 par step pmin pmax errpar
*
*	store fit result and error for A_bottom
*
	a_b    = par(3)
	a_berr = errpar(3)
	a_b_weight = [a_b]/([a_berr]*[a_berr])
	inv_a_berr = 1./([a_berr]*[a_berr])
	f_b    = par(4)
	f_berr = errpar(4)
	f_b_weight = [f_b]/([f_berr]*[f_berr])
	inv_f_berr = 1./([f_berr]*[f_berr])
*	r_b    = par(7)
*	r_berr = errpar(7)
*	r_b_weight = [r_b]/([r_berr]*[r_berr])
*	inv_r_berr = 1./([r_berr]*[r_berr])
*
*	if ( [opt] .eq. 'strexp' ) then
*		beta_b    = par(8)
*		beta_berr = errpar(8)
*		beta_b_weight = [beta_b]/([beta_berr]*[beta_berr])
*		inv_beta_berr = 1./([beta_berr]*[beta_berr])
*	endif
*
*
* plot result:

	pics kum
	set ysiz 25
	zone 2 3

	hi/pl [idl] e
	hi/pl [idt] e
	hi/pl [idr] e
	hi/pl [idb] e
*
*--------------------------------
*
*	calculate the weighted mean and the "inner" and "outer" standard deviation
*	fo asymmetry
*
	sum1 = [a_l_weight] + [a_t_weight] + [a_r_weight] + [a_b_weight]
	sum2 = [inv_a_lerr] + [inv_a_terr] + [inv_a_rerr] + [inv_a_berr]
	a_mean = [sum1]/[sum2]
	in_a_mean = $sigma(sqrt(1./[sum2]))  | inner standard deviation
*
*	now calculate the outer standard deviation
*
	dev_l = ([a_mean]-[a_l])*([a_mean]-[a_l])*[inv_a_lerr]
	dev_t = ([a_mean]-[a_t])*([a_mean]-[a_t])*[inv_a_terr]
	dev_r = ([a_mean]-[a_r])*([a_mean]-[a_r])*[inv_a_rerr]
	dev_b = ([a_mean]-[a_b])*([a_mean]-[a_b])*[inv_a_berr]
	sum3 = [dev_l] + [dev_t] + [dev_r] + [dev_b]
	out_a_mean = $sigma(sqrt(1./3.*[sum3]/[sum2]))
*
*----------------------------------
*	
*	mean of frequency
*
	sum1 = [f_l_weight] + [f_t_weight] + [f_r_weight] + [f_b_weight]
	sum2 = [inv_f_lerr] + [inv_f_terr] + [inv_f_rerr] + [inv_f_berr]
	f_mean = [sum1]/[sum2]
	in_f_mean = $sigma(sqrt(1./[sum2]))  | inner standard deviation
*
*	no calculate the outer standard deviation
*
	dev_l = ([f_mean]-[f_l])*([f_mean]-[f_l])*[inv_f_lerr]
	dev_t = ([f_mean]-[f_t])*([f_mean]-[f_t])*[inv_f_terr]
	dev_r = ([f_mean]-[f_r])*([f_mean]-[f_r])*[inv_f_rerr]
	dev_b = ([f_mean]-[f_b])*([f_mean]-[f_b])*[inv_f_berr]
	sum3 = [dev_l] + [dev_t] + [dev_r] + [dev_b]
	out_f_mean = $sigma(sqrt(1./3.*[sum3]/[sum2]))
*
*--------------------------------------
*
*	mean of relaxation
*
*	sum1 = [r_l_weight] + [r_t_weight] + [r_r_weight] + [r_b_weight]
*	sum2 = [inv_r_lerr] + [inv_r_terr] + [inv_r_rerr] + [inv_r_berr]
*	r_mean = [sum1]/[sum2]
*	in_r_mean = $sigma(sqrt(1./[sum2]))  | inner standard deviation
*
*	no calculate the outer standard deviation
*
*	dev_l = ([r_mean]-[r_l])*([r_mean]-[r_l])*[inv_r_lerr]
*	dev_t = ([r_mean]-[r_t])*([r_mean]-[r_t])*[inv_r_terr]
*	dev_r = ([r_mean]-[r_r])*([r_mean]-[r_r])*[inv_r_rerr]
*	dev_b = ([r_mean]-[r_b])*([r_mean]-[r_b])*[inv_r_berr]
*	sum3 = [dev_l] + [dev_t] + [dev_r] + [dev_b]
*	out_r_mean = $sigma(sqrt(1./3.*[sum3]/[sum2]))
*
*	mean of stretched exponent
*
*
*	if ( [opt] .eq. 'strexp' ) then
*		sum1 = [beta_l_weight] + [beta_t_weight] + [beta_r_weight] + [beta_b_weight]
*		sum2 = [inv_beta_lerr] + [inv_beta_terr] + [inv_beta_rerr] + [inv_beta_berr]
*		beta_mean = [sum1]/[sum2]
*		in_beta_mean = $sigma(sqrt(1./[sum2]))  | inner standard deviation
*
*		no calculate the outer standard deviation
*
*		dev_l = ([beta_mean]-[beta_l])*([beta_mean]-[beta_l])*[inv_beta_lerr]
*		dev_t = ([beta_mean]-[beta_t])*([beta_mean]-[beta_t])*[inv_beta_terr]
*		dev_r = ([beta_mean]-[beta_r])*([beta_mean]-[beta_r])*[inv_beta_rerr]
*		dev_b = ([beta_mean]-[beta_b])*([beta_mean]-[beta_b])*[inv_beta_berr]
*		sum3 = [dev_l] + [dev_t] + [dev_r] + [dev_b]
*		out_beta_mean = $sigma(sqrt(1./3.*[sum3]/[sum2]))
*
*	endif
*
*----------------------------------------
*
* output fit parameter definitions to text window:

	mess
	mess 'par(1): Background              [1]'
	mess 'par(2): normalization           [1]'
	mess 'par(3): asymmetry               [1]'
	mess 'par(4): frequency               [MHz]'
	mess 'par(5): phase                   [degree]'
	mess 'par(6): muon lifetime           [microSec]'
	mess 'par(7): relaxation              [1/microSec]'
	mess
	mess
	mess A = [a_mean]
	mess A_err_in  = [in_a_mean] 	
	mess A_err_out = [out_a_mean]
	mess 
	mess Freq. = [f_mean]
	mess F_err_in  = [in_f_mean] 	
	mess F_err_out = [out_f_mean]
	mess 
*	mess Relax. = [r_mean]
*	mess r_err_in  = [in_r_mean] 	
*	mess r_err_out = [out_r_mean]
*	mess 
*	if ( [opt] .eq. 'strexp' ) then
*		mess
*		mess Beta = [beta_mean]
*		mess Beta_err_in  = [in_beta_mean] 	
*		mess Beta_err_out = [out_beta_mean]
*		mess 
*	endif
	
* output to graphics window (quick and dirty)
	v/cre x(2) r 0 10
	v/cre y(2) r 0 10
	set xtic 0.0001
	set ytic 0.0001
	set vsiz 0.0001
	zone 1 3 3 s
	graph 2 x y w
*
	asym = $format([a_mean],F8.4)	
	in_err = $format([in_a_mean],F7.4)	
	out_err = $format([out_a_mean],F7.4)
*
	text 1 8 'A = '//[asym]//'+-'//[in_err]//'+-'//[out_err] 0.4
*
	freq = $format([f_mean],F8.5)
	in_err = $format([in_f_mean],F7.5)
	out_err = $format([out_f_mean],F7.5)
*	
	app = ' MHz'
	text 1 5 '[n] = '//[freq]//'+-'//[in_err]//'+-'//[out_err]//[app] 0.4


* restore previous settings:

	cd [current]
	zone
*	title ' '
	set xtic
	set ytic
	set vsiz
	set ysiz
	set xsiz