#!/usr/bin/env python

#########################################################
from os import makedirs, path				#
from os import remove, system, chdir			#
from numpy import array, transpose, dot         	#
from numpy.linalg import norm                   	#
from shutil import move, copy                   	#
from sys import stdout                         		# 
from math  import cos, sqrt, atan, fabs, sin, acos	#
########################################################
						#
first = 1					#
LAST =	1000					#
						#
#################################################

Reacted		= 0	# counters for raction probability
H		= 0	#
D		= 0	#
Scattered	= 0	#
Unclear		= 0	#
Running		= 0	# 
Trapped		= 0 	#
not_sub		= 0	#

Reacted_label	= [ ]	# tarj numbers
D_label 	= [ ]	#
H_label 	= [ ]	#
Unclear_label	= [ ]	#
Running_label	= [ ]	#
Trapped_label	= [ ]	#
################################################################################################################################################

last = 0

for i in range(first,(LAST+1)):		#cicle on the trajectories	

	dir = '%.4d' % i                # define workdir 
	if path.exists('%s' % dir):

		last += 1
		chdir(dir)                      # move in workdir	
	
		if not path.exists('outcome') and path.exists('POTCAR'):
			outcome = 'RUNNING'
			Running += 1 
			Running_label.append(dir)			

		elif not path.exists('outcome') and not path.exists('POTCAR'):
			outcome = 'not submitted'
			not_sub += 1 
		else:  	                       
			outcome = open('outcome').readline().split()[0]	

		if   outcome == 'REACTION'   : 
			Reacted  += 1
			Reacted_label.append(dir)
		 	diss = open('outcome').readline().split()[2]
			if diss == 'H':
				H += 1 
				H_label.append(dir)				
			else:
				D += 1
				D_label.append(dir)	

		elif outcome == 'UNCLEAR'  : 
			Unclear += 1
			Unclear_label.append(dir)
	
		elif outcome == 'SCATTERING':
			Scattered += 1
				
		elif outcome == 'TRAPPING':	
			Trapped += 1
			Trapped_label.append(dir)

		chdir('..')
	
	
### PROBABILITIES ################################################################################################################################

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

###################################################################################################################################################

Completed_trajs = last - Running - Unclear - not_sub

if Completed_trajs == 0 or Reacted == 0:  
	Completed_trajs = 999999
	Reacted		= 999999

Reacted_RP	= float(Reacted)			/ Completed_trajs 
React_Trap_RP	= float(Reacted+Trapped)		/ Completed_trajs 
H_RP		= float(H)				/ Completed_trajs
D_RP		= float(D)				/ Completed_trajs
F_CH            = float(H)                              / Reacted
R_confidence = 		confidence(Reacted_RP,	Completed_trajs)
RT_confidence =		confidence(React_Trap_RP,	Completed_trajs)
H_confidence =		confidence(H_RP,	Completed_trajs)
D_confidence =		confidence(D_RP,	Completed_trajs)
F_confidence =          confidence(F_CH,        Reacted)
### WRITE SUMMARY ################################################################################################################################

outSUM = open('Summary.dat', 'w')

outSUM.write('-----------------------------------\n')
outSUM.write('-- TRAJECTORIES ANALYSIS RESULTS --\n')
outSUM.write('-----------------------------------\n\n')

outSUM.write('Number of trajectories	= %.4d\n' % (last - first + 1))
outSUM.write('Of which, completed	= %.4d\n\n' % Completed_trajs)

outSUM.write('Reactive Events	= %d\n' % Reacted)
outSUM.write('\n')

outSUM.write('\nH dissociations	= %d\n' % H)
for j in range(len(H_label)): outSUM.write('%s    ' % H_label[j])
outSUM.write('\n')

outSUM.write('\nD dissociations	= %d\n' % D)
for j in range(len(D_label)): outSUM.write('%s    ' % D_label[j])

outSUM.write('\n\nTrapped trajs	= %d\n' % Trapped)
for j in range(len(Trapped_label)): outSUM.write('%s    ' % Trapped_label[j])

outSUM.write('\n\nScattering trajs= %d\n' % Scattered)


outSUM.write('-------------------------------------\n')

outSUM.write('Unclear trajs	= %d\n' % Unclear)
for j in range(len(Unclear_label)): outSUM.write('%s    ' % Unclear_label[j])
outSUM.write('\n\n')

outSUM.write('Running trajs	= %d\n' % Running)
for j in range(len(Running_label)): outSUM.write('%s	' % Running_label[j])
outSUM.write('\n')

outSUM.write('\n')
outSUM.write('-------------------------------------\n')

outSUM.write('		  S_0	\n')
outSUM.write('Total		%.3f +/- %.3f %%\n' 	% (Reacted_RP, R_confidence))
outSUM.write('Tot+Trap	%.3f +/- %.3f %%\n' 	% (React_Trap_RP, RT_confidence))
outSUM.write('\n-------------------------------------\n')
outSUM.write('H diss		%.3f +/- %.3f %%\n'     % (H_RP	, H_confidence))
outSUM.write('D diss		%.3f +/- %.3f %%\n'	% (D_RP	, D_confidence))
outSUM.write('\nFraction CH cleavage:  %.3f +/- %.3f %%\n'     % (F_CH, F_confidence))
outSUM.write('-------------------------------------\n')

outSUM.close()


print '\n'

Sum = open('Summary.dat').readlines()
for i in range(len(Sum)): print Sum[i].rstrip('\n')
print '\n'