#!/usr/bin/env python
import sys
from os import path, getcwd, popen
from os import chdir, system
import numpy as np

#folders = [ '0.94_v2j1', '1.18_v2j1', '1.29_v2j1', '1.55_v2j1', '1.80_v2j1', '2.12_v2j1', '2.56_v2j1' ]
#folders = [ '0.52_v1j1', '0.78_v1j1', '0.94_v1j1', '1.18_v1j1' ]
#folders = [ '0.28_ms_v1j1', '0.32_ms_v1j1',  '0.45_ms_v1j1',  '0.52_ms_v1j1',  '0.78_ms_v1j1',  '0.97_ms_v1j1',  '1.27_ms_v1j1' ]
#folders = [ '0.52_ms_v0j1', '0.78_ms_v0j1',  '0.97_ms_v0j1',  '1.27_ms_v0j1' ]
#folders = ['0.59_ms_v1j1', '0.64_ms_v1j1', '0.92_ms_v1j1', '0.94_ms_v1j1', '1.04_ms_v1j1']
folders = ['0.59_ms_v1j1', '0.64_ms_v1j1']

def confidence(P,T):
        'compute the 1 sigma confidece given a probability P and the total number of events T'
        return np.sqrt( P * (1-P) / (T-1) )

here=getcwd()
def collectvj(folder):
	chdir(folder)

	Ntraj = int( popen('ls | egrep "^[0-9].....$" | wc -l').read() )

	vj = []
	vj_weights = []
	here2=getcwd()				#workdir
	for i in range(1,(Ntraj+1)):
		traj= '%.6d' % i		# define working directory
		newdir=path.join(here2, traj)	# updating mypath
#		print( newdir )

		if not path.exists(newdir): 
			print( "folder", newdir, "not found. exiting..." )
#			continue
			quit()

		chdir(newdir)

		if path.exists('outcome'):
			outcome = open('outcome', 'r').readlines()
			if not 'SCATTERING\n' in outcome: continue

			lines = open( 'PostAnalysis.dat', 'r' ).readlines()
			if len(lines) < 2: continue
			if int( lines[33].split()[0] ) < 0 or int( lines[33].split()[1] ) < 0: continue
			vj.append( [ int( lines[32].split()[0] ), int( lines[32].split()[1] ), int( lines[33].split()[0] ), int( lines[33].split()[1] ) ] )
#			vj_weights.append( float( lines[37].split()[1] ) )	# 1GB nu, j, an, am
			vj_weights.append( float( lines[37].split()[1] ) )	# GB nu
#			vj_weights.append( float( lines[37].split()[1] ) * float( lines[37].split()[2] ) )	# GB nu, j
#			vj_weights.append( float( lines[37].split()[1] ) * float( lines[37].split()[2] ) * float( lines[37].split()[3] ) * float( lines[37].split()[4] ) )	# GB nu, j, an, am

	chdir(here)

	vj = np.array( vj )
	vj_weights = np.array( vj_weights )

	N = len(vj)
	Ngauss = sum( vj_weights[:] )

	vj_bin = np.bincount( vj[:,2] )

	print('------------------------------------------------------------\n')
	print( 'Histogram binning:' )
	print( 'Total weight: {0:f}'.format( N ) )

	print( '{:s}: Pvj->Pxj'.format( folder ) )
	print( vj_bin / float( N ) )
	if len( vj_bin ) > 2:
		print( '{:s}: Tvj->T1j +- SE, Tvj->T2j +- SE'.format( folder ) )
		T11 = vj_bin[1] / float( (vj_bin[1] + vj_bin[2]) )
		T11_SE = confidence( T11, N )
		T12 = vj_bin[2] / float( (vj_bin[1] + vj_bin[2]) )
		T12_SE = confidence( T12, N )
		print( '{:7.6f} {:7.6f} {:7.6f} {:7.6f}'.format( T11, T11_SE, T12, T12_SE ) )
	print( '{:s}: Tvj->T0j, Tvj->T1j'.format( folder ) )
	T00 = vj_bin[0] / float( (vj_bin[0] + vj_bin[1]) )
	T00_SE = confidence( T00, N )
	T01 = vj_bin[1] / float( (vj_bin[0] + vj_bin[1]) )
	T01_SE = confidence( T01, N )
	print( '{:7.6f} {:7.6f} {:7.6f} {:7.6f}'.format( T00, T00_SE, T01, T01_SE ) )

	vj_bin = np.bincount( vj[:,2], vj_weights[:] )

	print( '\n1Gaussian binning (nu, j, an, am):' )
	print( 'Total weight: {0:f}'.format( Ngauss ) )

	print( '{:s}: Pvj->Pxj'.format( folder ) )
	print( vj_bin / Ngauss )
	if len( vj_bin ) > 2:
		print( '{:s}: Tvj->T1j +- SE, Tvj->T2j +- SE'.format( folder ) )
		T11 = vj_bin[1] / float( (vj_bin[1] + vj_bin[2]) )
		T11_SE = confidence( T11, N )
		T12 = vj_bin[2] / float( (vj_bin[1] + vj_bin[2]) )
		T12_SE = confidence( T12, N )
		print( '{:7.6f} {:7.6f} {:7.6f} {:7.6f}'.format( T11, T11_SE, T12, T12_SE ) )
	print( '{:s}: Tvj->T0j, Tvj->T1j'.format( folder ) )
	T00 = vj_bin[0] / float( (vj_bin[0] + vj_bin[1]) )
	T00_SE = confidence( T00, N )
	T01 = vj_bin[1] / float( (vj_bin[0] + vj_bin[1]) )
	T01_SE = confidence( T01, N )
	print( '{:7.6f} {:7.6f} {:7.6f} {:7.6f}'.format( T00, T00_SE, T01, T01_SE ) )

for i in range(len(folders)):
	collectvj(folders[i])