source: fact/tools/pyscripts/simulation/absorb_photons.py

#!/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()