source: fact/tools/pyscripts/simulation/old_and_tests/check_corsi.py

#!/usr/bin/python -itt

import sys
import numpy as np

import rlcompleter
import readline
readline.parse_and_bind('tab: complete')

from ROOT import *

import readcorsika
from reflector import *


        
run = readcorsika.read_corsika_file(sys.argv[1])

h1 = TH2F("h_phi_theta","phi versus theta",180,0,180,360,-180,180)
huv = TH2F("huv"," u versus v",100,-1 ,1 ,100,-1,1)
huvm = TH2F("huvmean"," u-u_mean versus v-v_mean ",100,-0.3 ,0.3 ,100,-0.3,0.3)


event_counter = 0
for event in run.events:
    #print event_counter
    event_counter += 1
    if event.photons is None:
        continue
    for photon in event.photons:
        # photon is a 1D-nd.array with 10 elements
        # make Photon instance from this horrible array
        photon[1:3] -= np.array(event.info_dict['core_loc'])
        huv.Fill(photon[3],photon[4])
        photon = Photon(photon)
        d = photon.dir

        r = length(d)
        p = np.arctan2(d[1],d[0]) / np.pi * 180.
        t = np.arccos(d[2]/r) /np.pi * 180.

        h1.Fill( t, p )
        
        
        
c1 = TCanvas("c1","title",0,0,500,500)
h1.Draw("colz")
c1.Update()

c2 = TCanvas("c2","title",500,0,500,500)
huv.Draw("colz")
c2.Update()

print "u mean:",huv.GetMean(1)
print "v mean",huv.GetMean(2)
umean = huv.GetMean(1)
vmean = huv.GetMean(2)

for event in run.events:
    #print event_counter
    event_counter += 1
    if event.photons is None:
        continue
    for photon in event.photons:
        
        huvm.Fill(photon[3]-umean,photon[4]-vmean)
        
c3 = TCanvas("c3","title",0,500,500,500)
huvm.Draw("colz")
c3.Update()