#!/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 reaction probability
OH		= 0	#
CH		= 0	#
CO		= 0	#
Scattered	= 0	#
Unclear		= 0	#
Running		= 0	# 
Trapped		= 0 	#
not_sub		= 0	#

Reacted_label	= [ ]	# tarj numbers
OH_label 	= [ ]	#
CH_label 	= [ ]	#
CO_label 	= [ ]	#
Unclear_label	= [ ]	#
Running_label	= [ ]	#
Trapped_label	= [ ]	#
################################################################################################################################################

last = 0

for i in range(first,(LAST+1)):		#cycle 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 == 'OH':
				OH += 1 
				OH_label.append(dir)				
			elif diss == 'CH1' or diss == 'CH2' or diss == 'CH3':
				CH += 1
				CH_label.append(dir)
			else:
				CO += 1
				CO_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 
OH_RP		= float(OH)				/ Completed_trajs
CH_RP		= float(CH)				/ Completed_trajs
CO_RP		= float(CO)				/ Completed_trajs
F_OH            = float(OH)                             / Reacted
F_CH            = float(CH)                             / Reacted
F_CO            = float(CO)                             / Reacted
R_confidence = 		confidence(Reacted_RP,	Completed_trajs)
RT_confidence =		confidence(React_Trap_RP,	Completed_trajs)
OH_confidence =		confidence(OH_RP,	Completed_trajs)
CH_confidence =		confidence(CH_RP,	Completed_trajs)
CO_confidence =		confidence(CO_RP,	Completed_trajs)
F_OH_confidence =       confidence(F_OH,        Reacted)
F_CH_confidence =       confidence(F_CH,        Reacted)
F_CO_confidence =       confidence(F_CO,        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('\nOH dissociations	= %d\n' % OH)
for j in range(len(OH_label)): outSUM.write('%s    ' % OH_label[j])
outSUM.write('\n')

outSUM.write('\nCH dissociations	= %d\n' % CH)
for j in range(len(CH_label)): outSUM.write('%s    ' % CH_label[j])
outSUM.write('\n')

outSUM.write('\nCO dissociations	= %d\n' % CO)
for j in range(len(CO_label)): outSUM.write('%s    ' % CO_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('OH diss		%.3f +/- %.3f %%\n'     % (OH_RP	, OH_confidence))
outSUM.write('CH diss		%.3f +/- %.3f %%\n'	% (CH_RP	, CH_confidence))
outSUM.write('CO diss		%.3f +/- %.3f %%\n'     % (CO_RP        , CO_confidence))
outSUM.write('\nFraction OH cleavage:  %.3f +/- %.3f %%\n'     % (F_OH, F_OH_confidence))
outSUM.write('\nFraction CH cleavage:  %.3f +/- %.3f %%\n'     % (F_CH, F_CH_confidence))
outSUM.write('\nFraction CO cleavage:  %.3f +/- %.3f %%\n'     % (F_CO, F_CO_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'