-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathSimulatedEvents.jl
250 lines (226 loc) · 10.2 KB
/
SimulatedEvents.jl
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
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
using ROOTFramework, Cxx
using Clustering
using MJDSigGen: fieldgen, signal_calc_init, get_signal!, drift_path, read_fields, outside_detector, nearest_field_grid_index, cyl2cart, cart2cyl
type Event
hitpositions
hitenergies
hitpulses
hitepaths
hithpaths
hitclusters
clusterpositions
clustersigmas
clusterenergies
clusterpulses
clusterepaths
clusterhpaths
totenergy
sumpulse
end
#function to read in converted MaGe events and generate pulses from them
function create_events(setup::MJDSigGen.Struct_MJD_Siggen_Setup, root_file_name, root_tree_name; min_energy = 0, max_energy = 9999, save_hits = false, save_clusters = false, cluster_size = -1f0)
#input is a converted MaGe root file that has a hittree
#min_energy = 600 is useful to save memory and storage space
#save_hits = true saves all hits: positions, energies, pulses, drift paths
#save_clusters = true saves all clusters: positions, sizes, energies, pulses, drift paths
#cluster_size <= 0 or Nhits <= 3 the sum of hit pulses will be saved, otherwise the sum of cluster pulses will be saved
#clustering only works for more than 3 hits :S
if save_clusters && cluster_size <= 0
error("Error: cannot save clusters with size <= 0")
end
#root file input
bindings = TTreeBindings()
Nhits = bindings[:Nhits] = Ref(zero(Int32))
TotEnergy = bindings[:TotEnergy] = Ref(zero(Float32))
HitEnergy = bindings[:HitEnergy] = ROOTFramework.CxxObjWithPtrRef(icxx"std::vector<float>();")
HitX = bindings[:HitX] = ROOTFramework.CxxObjWithPtrRef(icxx"std::vector<float>();")
HitY = bindings[:HitY] = ROOTFramework.CxxObjWithPtrRef(icxx"std::vector<float>();")
HitZ = bindings[:HitZ] = ROOTFramework.CxxObjWithPtrRef(icxx"std::vector<float>();")
#initialize empty output
allevents = Array{Event}(0)
#loop over root file
open(TChainInput, bindings, root_tree_name, root_file_name) do input
i = 0
for _ in input
i += 1
if TotEnergy.x < min_energy
continue
end
if TotEnergy.x > max_energy
continue
end
#initializations
eventpulse = zero(Float32)
goodevent = true #workaround for events with hits at the contact or groove
setup.charge_cloud_size = 0.0
setup.use_diffusion = 1
if cluster_size > 0 || save_hits || save_clusters
positions = zeros(Float32, 3, Nhits.x)
energies = zeros(Float32, Nhits.x)
if save_hits
pulses = Array{Array{Float32}}(0)
epaths = Array{Array{Float32,2}}(0)
hpaths = Array{Array{Float32,2}}(0)
end
if save_clusters
clpositions = Array{Array{Float32}}(0)
clsigmas = Array{Float32}(0)
clenergies = Array{Float32}(0)
clpulses = Array{Array{Float32}}(0)
clepaths = Array{Array{Float32,2}}(0)
clhpaths = Array{Array{Float32,2}}(0)
end
end
#loop over hits
for j in 0:(Nhits.x-1)
if cluster_size > 0 || save_hits || save_clusters
positions[1,j+1], positions[2,j+1], positions[3,j+1] = HitX.x[j], HitY.x[j], HitZ.x[j]
energies[j+1] = HitEnergy.x[j]
end
if cluster_size <= 0 || Nhits.x <= 3 || save_hits
p1 = (HitX.x[j], HitY.x[j], HitZ.x[j])
ind = nearest_field_grid_index(setup, p1)
setup.energy = HitEnergy.x[j]
if ind[1] != :outside
if ind[1] == :interpol
pulse1 = get_signal!(setup, p1)
elseif ind[1] == :extrapol
#println("event $i point need extrapolation $p1")
phi = cart2cyl(p1[1],p1[2],p1[3])[2]
p1 = cyl2cart(ind[2]*setup.xtal_grid, phi, ind[3]*setup.xtal_grid)
pulse1 = get_signal!(setup, p1)
else
println("this should not happen")
end
if cluster_size <= 0 || Nhits.x <= 3
eventpulse += HitEnergy.x[j] * pulse1
end
if save_hits
push!(pulses, HitEnergy.x[j] * pulse1)
push!(epaths, drift_path(setup, :e))
push!(hpaths, drift_path(setup, :h))
end
else
println("event $i skipped because point outside $p1")
goodevent = false
end
end
end
#do clustering
if cluster_size > 0 && Nhits.x > 3
clusters = dbscan(positions, cluster_size, min_neighbors=1, min_cluster_size=1, leafsize=20)
for j in 1:length(clusters)
inds = vcat(clusters[j].core_indices, clusters[j].boundary_indices)
Ecl = sum(energies[inds])
p2 = [sum(energies[inds].*positions[1,inds]); sum(energies[inds].*positions[2,inds]); sum(energies[inds].*positions[3,inds])]./Ecl
pulse2 = zero(Float32)
check = nearest_field_grid_index(setup, tuple(p2...))
if check[1] != :outside
clsize = 0
for k=1:length(inds)
clsize += energies[inds[k]] * (norm(positions[:,inds[k]]-p2)^2)
end
clsize = sqrt(clsize/Ecl)
setup.charge_cloud_size = clsize
setup.energy = Ecl
if check[1] == :interpol
pulse2 = get_signal!(setup, tuple(p2...))
elseif check[1] == :extrapol
#println("event $i cluster need extrapolation $p2")
phi = cart2cyl(p2[1],p2[2],p2[3])[2]
p2[1],p2[2],p2[3] = cyl2cart(check[2]*setup.xtal_grid, phi, check[3]*setup.xtal_grid)
checkagain = nearest_field_grid_index(setup, tuple(p2...))
if checkagain[1] == :interpol
pulse2 = get_signal!(setup, tuple(p2...))
else
println("event $i skipped because extrapolated cluster outside $p2")
goodevent = false
end
else
println("this should not happen")
end
if save_clusters && goodevent
push!(clenergies, Ecl)
push!(clpositions, p2)
push!(clsigmas, clsize)
push!(clpulses, Ecl * pulse2)
push!(clepaths, drift_path(setup, :e))
push!(clhpaths, drift_path(setup, :h))
end
eventpulse += Ecl * pulse2
else
println("event $i skipped because cluster outside $p2")
goodevent = false
end
end
end
#save everything
if goodevent
event1 = Event(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, TotEnergy.x, eventpulse)
if save_hits
event1.hitpositions = positions
event1.hitenergies = energies
event1.hitpulses = pulses
event1.hitepaths = epaths
event1.hithpaths = hpaths
if cluster_size > 0 && Nhits.x > 3
event1.hitclusters = clusters
end
end
if save_clusters
event1.clusterpositions = clpositions
event1.clustersigmas = clsigmas
event1.clusterenergies = clenergies
event1.clusterpulses = clpulses
event1.clusterepaths = clepaths
event1.clusterhpaths = clhpaths
end
push!(allevents, event1)
end
end #loop over tree
end #loop over root file
allevents
end
function plot_crystal_3d(setup::MJDSigGen.Struct_MJD_Siggen_Setup)
R = setup.xtal_radius
L = setup.xtal_length
LT = setup.outer_taper_length
WT = setup.outer_taper_width
cyl_x = cos.(linspace(0,2π,40))
cyl_y = sin.(linspace(0,2π,40))
cyl_z = zeros(cyl_x)
myfig = plot(cyl_x*R, cyl_y*R, cyl_z,
label = "", l = (:grey, :dot),
xaxis = ("x (mm)", (-R, R)), yaxis = ("y (mm)", (-R, R)), zaxis = ("z (mm)", (0, L)),
size = (600,600),
)
for z in linspace(0,L-LT,5)
plot!(cyl_x*R, cyl_y*R, cyl_z+z, label = "", linecolor = :grey, line = :dot)
end
for z in linspace(L-LT,L,5)
r = -z*WT/LT + WT/LT*L + R-WT
plot!(cyl_x*r, cyl_y*r, cyl_z+z, label = "", linecolor = :grey, line = :dot)
end
for r in linspace(5,R,5)
plot!(cyl_x*r, cyl_y*r, cyl_z, label = "", linecolor = :grey, line = :dot)
end
for r in linspace(5,R-WT,5)
plot!(cyl_x*r, cyl_y*r, cyl_z+L, label = "", linecolor = :grey, line = :dot)
end
myfig
end
#function to write simple root tree with energy and amplitude values
function write_aoe_tree(allevents::Array{Event}, root_file_name)
bindings = TTreeBindings()
Energy = bindings[:Energy] = Ref(zero(Float64))
Amplitude = bindings[:Amplitude] = Ref(zero(Float64))
open(TFile, root_file_name, "recreate") do tfile
ttree = create_ttree!(tfile, "simulation", "Simulated events")
output = TTreeOutput(ttree, bindings)
for i in 1:length(allevents)
Energy.x = allevents[i].totenergy
Amplitude.x = maximum(vcat([0],diff(allevents[i].sumpulse)))
push!(output)
end
end
end