#!/usr/bin/env python
import sys, commands
from os import path, getcwd
from os import chdir, system
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

fig, ax = plt.subplots()

first = 1
last  = 500

r_ts = 1.61

Zcom_ts = []
here=getcwd()				#workdir
for i in range(first,(last+1)):
	traj= '%.4d' % i		# define working directory
	newdir=path.join(here, traj)	# updating mypath

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

	chdir(newdir)

	if path.exists('analysis_com.dat') and path.exists('analysis_r.dat') and path.exists('outcome'):
		outcome  = open('outcome', 'r').readline()
	        if "SCATTERING" in outcome or "TRAPPING" in outcome:
			continue
		elif "REACTION" in outcome:
			com_dat  = open('analysis_com.dat', 'r').readlines()
			r_dat = open('analysis_r.dat', 'r').readlines()
			breakloop = False
			for j in range(2, len(r_dat)):
				for k in range(1,5):
					if float(r_dat[j].split()[k]) > r_ts:
						Zcom_ts.append( float(com_dat[j].split()[6]) )
						breakloop = True
						break
				if breakloop == True: break

n, bins, patches = plt.hist(Zcom_ts, 10, normed=True)
print np.mean(Zcom_ts)
print np.std(Zcom_ts) / (len(Zcom_ts)-1)**0.5

plt.show()