#!/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 = 0					#
LAST =	100000					#
						#
#################################################

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

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

last = 0

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

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

		last += 1
		chdir(dir)                      # move in workdir	

		if path.exists('BONDveryLONG'):
			Verylong += 1
			Verylong_label.append(dir)
	
		if not path.exists('outcome'):
			outcome = 'RUNNING'
			Running += 1 
			Running_label.append(dir)
		else:  	                       
			outcome = open('outcome').readline().split()[0]	

			if 'REACTION' in outcome: 
				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 'SCATTERING' in outcome:
				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

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	= %.6d\n' % (last - first))
outSUM.write('Of which, completed	= %.6d\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')

outSUM.write('\nTrajs with very long bonds	= %d\n' % Trapped)
for j in range(len(Verylong_label)): outSUM.write('%s    ' % Verylong_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('\n')
outSUM.write('-------------------------------------\n')

outSUM.write('		  S_0	\n')
outSUM.write('Total		%.6f +/- %.6f %%\n' 	% (Reacted_RP, R_confidence))
outSUM.write('Tot+Trap	%.6f +/- %.6f %%\n' 	% (React_Trap_RP, RT_confidence))
outSUM.write('\n-------------------------------------\n')
outSUM.write('H diss		%.6f +/- %.6f %%\n'     % (H_RP	, H_confidence))
outSUM.write('D diss		%.6f +/- %.6f %%\n'	% (D_RP	, D_confidence))
outSUM.write('\nFraction CH cleavage:  %.6f +/- %.6f %%\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' )