Analysis.py
4.17 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
#!/bin/python
from sys import argv
from numpy import append, savetxt, shape, array, newaxis, zeros, arange
from Modules.Read import ReadEnergy, ReadTime, ReadExtraFile, ReadProfile
from Modules.Spectrum import spectrum
from Modules.Angle import radial, angle_vs_energy
from Modules.Timing import timing
def error(error_type):
print error_type
exit()
if shape(argv)[0] < 2:
error("not enough arguments (at least 1)")
folder = "Results/"
for fileId in argv[1:]:
# read files
time = ReadTime(fileId)
energy = ReadEnergy(fileId)
weightini, generation, theta_arrival, Esource = ReadExtraFile(fileId,[2,3,4,5])
nbPhotonsEmitted=ReadProfile(fileId,[3])
Gen_contrib = []
Gen_cont = zeros((int(max(generation))+1))
Source = []
Spectrum = []
Radial = []
Angle_Energy = []
Timing = []
for powerlaw_index in [1,1.5,2,2.5]:
# apply source spectrum
weight = weightini*(Esource)**(1-powerlaw_index)
#=============================================================================#
# GENERATION CONTRIBUTION
#=============================================================================#
NbTotEvents = sum(weight)
if Gen_contrib==[]:
Gen_contrib = arange(0,max(generation)+1,1)[:,newaxis]
for gen in list(set(generation)):
Gen_cont[int(gen)] = sum(weight[generation==gen])/NbTotEvents *100
Gen_contrib = append(Gen_contrib,Gen_cont[:,newaxis],axis=1)
#=============================================================================#
# SPECTRUM (SOURCE AND MEASURED)
#=============================================================================#
nbBins = 100
# draw source spectrum
Es=array(list(set(Esource)))
Ws= (Es)**(1-powerlaw_index)
Es,Fs = spectrum(Es,Ws,nbBins)
Es=Es[:,newaxis]
Fs=Fs[:,newaxis]
if Source==[]:
Source = Es
Source = append(Source,Fs,axis=1)
else:
Source = append(Source,Fs,axis=1)
# primary gamma-rays contribution
cond = (generation==0)
ener,flux_0 = spectrum(energy[cond],weight[cond],nbBins)
# full spectrum
ener,flux = spectrum(energy,weight,nbBins)
ener=ener[:,newaxis]
flux=flux[:,newaxis]
flux_0=flux_0[:,newaxis]
if Spectrum==[]:
Spectrum = ener
Spectrum = append(Spectrum,flux,axis=1)
Spectrum = append(Spectrum,flux_0,axis=1)
else:
Spectrum = append(Spectrum,flux,axis=1)
Spectrum = append(Spectrum,flux_0,axis=1)
#=============================================================================#
# RADIAL DISTRIBUTION AND ANGLE VERSUS ENERGY
#=============================================================================#
nbBins = 100
ener,angle = angle_vs_energy(theta_arrival,energy,weight,nbPhotonsEmitted)
ener=ener[:,newaxis]
angle=angle[:,newaxis]
if Angle_Energy==[]:
Angle_Energy = ener
Angle_Energy = append(Angle_Energy,angle,axis=1)
else:
Angle_Energy = append(Angle_Energy,angle,axis=1)
theta,dndtheta2 = radial(theta_arrival,weight,nbPhotonsEmitted)
theta=theta[:,newaxis]
dndtheta2=dndtheta2[:,newaxis]
if Radial==[]:
Radial = theta
Radial = append(Radial,dndtheta2,axis=1)
else:
Radial = append(Radial,dndtheta2,axis=1)
#=============================================================================#
# TIME DISTRIBUTION AND TIME DELAY VERSUS ANGLE
#=============================================================================#
nbBins = 100
delta_t,dNdt = timing(time,weight,nbBins)
delta_t=delta_t[:,newaxis]
dNdt=dNdt[:,newaxis]
if Timing==[]:
Timing = delta_t
Timing = append(Timing,dNdt,axis=1)
else:
Timing = append(Timing,dNdt,axis=1)
savetxt(folder+fileId+"/Generation.txt",Gen_contrib)
savetxt(folder+fileId+"/Spectrum.txt",Spectrum)
savetxt(folder+fileId+"/Source_spectrum.txt",Source)
savetxt(folder+fileId+"/Angle_vs_Energy.txt",Angle_Energy)
savetxt(folder+fileId+"/Radial_distribution.txt",Radial)
savetxt(folder+fileId+"/Timing.txt",Timing)