Index: fact/tools/pyscripts/simulation/absorb_photons.py
===================================================================
--- fact/tools/pyscripts/simulation/absorb_photons.py	(revision 14820)
+++ fact/tools/pyscripts/simulation/absorb_photons.py	(revision 14820)
@@ -0,0 +1,67 @@
+#!/usr/bin/python -itt
+
+import sys
+import numpy as np
+import rlcompleter
+import readline
+readline.parse_and_bind('tab: complete')
+
+from readcorsika import read_corsika_file
+import fact_spline_interpolations as fact_si
+from Turnable import *
+
+absorbtion_coeff = fact_si.interpolate_from_files( 
+                fact_si.Wavelength_dependency_description_files )
+
+angular_acceptance = fact_si.interpolate_from_files( 
+                fact_si.angular_acceptance_files)
+
+
+run = read_corsika_file(sys.argv[1])
+
+print "photons = ", run.table.shape[0]
+
+wavelength = run.table[:,9]
+efficiencies = absorbtion_coeff(wavelength)
+# get uniform randomvalues between 0.0 and 1.0
+rndm = np.random.rand(wavelength.shape[0])
+
+# if the randomvalue is larger than the efficiency, the photon is absorbed
+# if the randomvalue is smaller, the photon survives
+run.not_absorbed = run.table[rndm < efficiencies,:]
+
+print "photons not absorbed = ", run.not_absorbed.shape[0], 
+print float(run.not_absorbed.shape[0])/run.table.shape[0]*100, "%"
+
+distance = np.sqrt((run.not_absorbed[:,0:3]**2).sum(axis=1))
+run.inside = run.table[distance<210,:]
+
+print "photons inside = ", run.inside.shape[0],
+print float(run.inside.shape[0])/run.not_absorbed.shape[0]*100, "%"
+
+tm = run.inside[:,6].mean()
+pm = run.inside[:,7].mean()
+print "theta mean = ", tm 
+print "  phi mean = ", pm
+# now I can turn the photons such that the photons come from directly above
+
+# lets turn the world around the z-axis, so that phi mean becomes zero.
+Rphi = make_rotation_matrix( [0,0,1], pm)
+
+new_cosines = np.dot( run.inside[:,3:6],Rphi)
+new_theta = np.arccos(new_cosines[:,2])
+new_phi = np.arctan2(new_cosines[:,1], new_cosines[:,0])
+
+print "new theta mean", new_theta.mean()
+print "new phi mean", new_phi.mean()
+
+# now I think we need to turn around the y-axis
+Rtheta = make_rotation_matrix( [0,1,0], tm)
+
+new_cosines = np.dot( new_cosines,Rtheta)
+new_theta = np.arccos(new_cosines[:,2])
+new_phi = np.arctan2(new_cosines[:,1], new_cosines[:,0])
+
+print "new theta mean", new_theta.mean()
+print "new phi mean", new_phi.mean()
+