#!/bin/bash
echo "test14: Eloss and multiple scattering in LH2 (6 sec)"

cat >test14.in <<-!
	param histoFile=""
	physics QGSP
	beam gaussian meanMomentum=200 particle=mu+ nEvents=10000
	trackcuts keep=mu+
	tubs Absorber length=1000 outerRadius=500 material=LH2
	place Absorber z=500
	virtualdetector AF format=BLTrackFile radius=500 length=1 color=1,1,1 file=test14.txt
	place AF z=1000.5
!

cat >test14.check <<-\!
	#!/usr/bin/tclsh
	# reads stdin, computing mean and sigma for X, Y, Ptot
	# exits successfully if both mean and sigma are within
	# tolerance of the good values (whcih came from a good run).
	# In any case, the means and sigmas are written to stdout.
	
	# columns of BLTrackfile (starts with 0)
	set c_x 0
	set c_y 1
	set c_px 3
	set c_py 4
	set c_pz 5
	set c_pid 7
	# tolerance is in percent
	set tolerance 10.0
	
	set sum_x 0.0
	set sum_y 0.0
	set sum_ptot 0.0
	set sumsq_x 0.0
	set sumsq_y 0.0
	set sumsq_ptot 0.0

	# read stdin, totaling sum and sumsq of ptot
	set input [split [read stdin] "\n"]
	set n 0
	foreach line $input {
		if {[string index $line 0] == "#"} { continue }
		set fields [split $line " "]
		set x [lindex $fields $c_x]
		set y [lindex $fields $c_y]
		set px [lindex $fields $c_px]
		set py [lindex $fields $c_py]
		set pz [lindex $fields $c_pz]
		if {"$x" == "" || "$y" == "" || "$px" == "" || "$py" == "" ||
		    "$py" == ""} {
			continue
		}
		set n [expr $n+1]
		set ptot [expr sqrt($px*$px+$py*$py+$pz*$pz)]
		set sum_x [expr $sum_x+$x]
		set sumsq_x [expr $sumsq_x+$x*$x]
		set sum_y [expr $sum_y+$y]
		set sumsq_y [expr $sumsq_y+$y*$y]
		set sum_ptot [expr $sum_ptot+$ptot]
		set sumsq_ptot [expr $sumsq_ptot+$ptot*$ptot]
	}
	if {$n<=1} { exit 99 }
	
	# convert sum and sumsq to mean and sigma
	set tcl_precision 3
	set mean_x [expr $sum_x/$n]
	set sigma_x [expr sqrt(abs($sumsq_x/$n - $mean_x*$mean_x))]
	set mean_y [expr $sum_y/$n]
	set sigma_y [expr sqrt(abs($sumsq_y/$n - $mean_y*$mean_y))]
	set mean_ptot [expr $sum_ptot/$n]
	set sigma_ptot [expr sqrt(abs($sumsq_ptot/$n - $mean_ptot*$mean_ptot))]
	puts "x: $mean_x+-$sigma_x y: $mean_y+-$sigma_y Ptot: $mean_ptot+-$sigma_ptot"
	
	# get good values
	#x: 0.139+-13.4 y: -0.0392+-14.6 Ptot: 163.0+-8.1
	set good_sigma_x 13.0
	set good_sigma_y 14.0
	set good_mean_ptot 163.0
	set good_sigma_ptot 6
	
	# verify current values are within tolerance of good values
	set tol [expr $tolerance/100.0]
	if {[expr abs($mean_x)>$good_sigma_x*$tol] ||
	    [expr abs($sigma_x-$good_sigma_x)>$good_sigma_x*$tol]} \
		{ puts "BAD X"; exit 1}
	if {[expr abs($mean_y)>$good_sigma_y*$tol] ||
	    [expr abs($sigma_y-$good_sigma_y)>$good_sigma_y*$tol]} \
		{ puts "BAD Y"; exit 1}
	if {[expr abs($mean_ptot-$good_mean_ptot)>$good_sigma_ptot*$tol] ||
	    [expr abs($sigma_ptot-$good_sigma_ptot)>$good_sigma_ptot*$tol]} \
		{ puts "BAD Ptot"; exit 1}
	
	exit 0
!
chmod +x test14.check

if ! g4bl test14.in >test14.out
then
	echo "****** test14: g4bl FAILED! -- test14.* files left intact"
	exit 1
fi

if ! ./test14.check <test14.txt >>test14.out
then
	echo "****** test14: mean & sigma test FAILED -- test14.* files left intact"
	exit 1
fi

# clean up and exit
rm test14.*
exit 0