135 def generator(nuslice_tree, rootfile, pset):
136 drift_distance = pset.DriftDistance
138 xbin_width = drift_distance/x_bins
139 half_bin_width = xbin_width/2.
140 y_bias_slope = pset.YBiasSlope
141 z_bias_slope = pset.ZBiasSlope
143 xvals = np.arange(half_bin_width, drift_distance, xbin_width)
144 xerrs = np.array([half_bin_width] * len(xvals))
145 dist_to_anode_bins = x_bins
146 dist_to_anode_low = 0.
147 dist_to_anode_up = drift_distance
148 if detector ==
"sbnd":
151 elif detector ==
"icarus":
152 if pset.Cryostat == 0:
155 if pset.Cryostat == 1:
159 profile_bins = x_bins
162 dy_spreads = [
None] * x_bins
163 dy_means = [
None] * x_bins
164 dy_h2 = TH2D(
"dy_h2",
"#Delta y",
165 dist_to_anode_bins, dist_to_anode_low, dist_to_anode_up,
166 pset.dy_bins, pset.dy_low, pset.dy_up)
167 dy_h2.GetXaxis().SetTitle(
"distance from anode (cm)")
168 dy_h2.GetYaxis().SetTitle(
"y_flash - y_TPC (cm)")
169 dy_prof = TProfile(
"dy_prof",
"Profile of dy_spreads in #Delta y",
170 profile_bins, dist_to_anode_low, dist_to_anode_up,
171 pset.dy_low*2, pset.dy_up*2, profile_option)
172 dy_prof.GetXaxis().SetTitle(
"distance from anode (cm)")
173 dy_prof.GetYaxis().SetTitle(
"y_flash - y_TPC (cm)")
174 dy_h1 = TH1D(
"dy_h1",
"",
175 profile_bins, dist_to_anode_low, dist_to_anode_up)
176 dy_h1.GetXaxis().SetTitle(
"distance from anode (cm)")
177 dy_h1.GetYaxis().SetTitle(
"y_flash - y_TPC (cm)")
179 dz_spreads = [
None] * x_bins
180 dz_means = [
None] * x_bins
181 dz_h2 = TH2D(
"dz_h2",
"#Delta z",
182 dist_to_anode_bins, dist_to_anode_low, dist_to_anode_up,
183 pset.dz_bins, pset.dz_low, pset.dz_up)
184 dz_h2.GetXaxis().SetTitle(
"distance from anode (cm)")
185 dz_h2.GetYaxis().SetTitle(
"z_flash - z_TPC (cm)")
186 dz_prof = TProfile(
"dz_prof",
"Profile of dz_spreads in #Delta z",
187 profile_bins, dist_to_anode_low, dist_to_anode_up,
188 pset.dz_low*2.5, pset.dz_up*2.5, profile_option)
189 dz_prof.GetXaxis().SetTitle(
"distance from anode (cm)")
190 dz_prof.GetYaxis().SetTitle(
"z_flash - z_TPC (cm)")
191 dz_h1 = TH1D(
"dz_h1",
"",
192 profile_bins, dist_to_anode_low, dist_to_anode_up)
193 dz_h1.GetXaxis().SetTitle(
"distance from anode (cm)")
194 dz_h1.GetYaxis().SetTitle(
"z_flash - z_TPC (cm)")
196 rr_spreads = [
None] * x_bins
197 rr_means = [
None] * x_bins
198 rr_h2 = TH2D(
"rr_h2",
"PE Spread",
199 dist_to_anode_bins, dist_to_anode_low, dist_to_anode_up,
200 pset.rr_bins, pset.rr_low, pset.rr_up)
201 rr_h2.GetXaxis().SetTitle(
"distance from anode (cm)")
202 rr_h2.GetYaxis().SetTitle(
"spread (cm)")
203 rr_prof = TProfile(
"rr_prof",
"Profile of PE Spread",
204 profile_bins, dist_to_anode_low, dist_to_anode_up,
205 pset.rr_low, pset.rr_up, profile_option)
206 rr_prof.GetXaxis().SetTitle(
"distance from anode (cm)")
207 rr_prof.GetYaxis().SetTitle(
"spread (cm)")
208 rr_h1 = TH1D(
"rr_h1",
"",
209 profile_bins, dist_to_anode_low, dist_to_anode_up)
210 rr_h1.GetXaxis().SetTitle(
"distance from anode (cm)")
211 rr_h1.GetYaxis().SetTitle(
"spread (cm)")
213 ratio_spreads = [
None] * x_bins
214 ratio_means = [
None] * x_bins
215 ratio_h2 = TH2D(
"ratio_h2",
"PE Ratio",
216 dist_to_anode_bins, dist_to_anode_low, dist_to_anode_up,
217 pset.ratio_bins, pset.ratio_low, pset.ratio_up)
218 ratio_h2.GetXaxis().SetTitle(
"distance from anode (cm)")
219 ratio_h2.GetYaxis().SetTitle(
"ratio")
220 ratio_prof = TProfile(
"ratio_prof",
"Profile of PE Ratio",
221 profile_bins, dist_to_anode_low, dist_to_anode_up,
222 pset.ratio_low, pset.ratio_up, profile_option)
223 ratio_prof.GetXaxis().SetTitle(
"distance from anode (cm)")
224 ratio_prof.GetYaxis().SetTitle(
"ratio")
225 ratio_h1 = TH1D(
"ratio_h1",
"",
226 profile_bins, dist_to_anode_low, dist_to_anode_up)
227 ratio_h1.GetXaxis().SetTitle(
"distance from anode (cm)")
228 ratio_h1.GetYaxis().SetTitle(
"ratio")
230 unfolded_score_scatter = TH2D(
"unfolded_score_scatter",
"Scatter plot of match scores",
231 2*dist_to_anode_bins, x_gl_low, x_gl_up,
232 pset.score_hist_bins, pset.score_hist_low, pset.score_hist_up*(3./5.))
233 unfolded_score_scatter.GetXaxis().SetTitle(
"X global (cm)")
234 unfolded_score_scatter.GetYaxis().SetTitle(
"match score (arbitrary)")
235 oldunfolded_score_scatter = TH2D(
"oldunfolded_score_scatter",
"Scatter plot of match scores, old metrics",
236 2*dist_to_anode_bins, x_gl_low, x_gl_up,
237 pset.score_hist_bins, pset.score_hist_low, pset.score_hist_up*(3./5.))
238 oldunfolded_score_scatter.GetXaxis().SetTitle(
"X global (cm)")
239 oldunfolded_score_scatter.GetYaxis().SetTitle(
"match score (arbitrary)")
240 match_score_scatter = TH2D(
"match_score_scatter",
"Scatter plot of match scores",
241 dist_to_anode_bins, dist_to_anode_low, dist_to_anode_up,
242 pset.score_hist_bins, pset.score_hist_low, pset.score_hist_up*(3./5.))
243 match_score_scatter.GetXaxis().SetTitle(
"distance from anode (cm)")
244 match_score_scatter.GetYaxis().SetTitle(
"match score (arbitrary)")
245 match_score_h1 = TH1D(
"match_score",
"Match Score",
246 pset.score_hist_bins, pset.score_hist_low, pset.score_hist_up)
247 match_score_h1.GetXaxis().SetTitle(
"match score (arbitrary)")
250 for e
in nuslice_tree:
251 if e.slices > e.true_nus:
continue
253 rr_h2.Fill(qX, e.flash_rr)
254 rr_prof.Fill(qX, e.flash_rr)
255 ratio_h2.Fill(qX, e.flash_ratio)
256 ratio_prof.Fill(qX, e.flash_ratio)
258 for e
in nuslice_tree:
259 if e.slices > e.true_nus:
continue
262 corr_flash_y = e.flash_yb -
y_bias(e.y_skew, hypo_x, y_bias_slope)
263 corr_flash_z = e.flash_zb -
z_bias(e.z_skew, hypo_x, z_bias_slope)
265 dy_h2.Fill(qX, corr_flash_y - e.charge_y)
266 dy_prof.Fill(qX, corr_flash_y - e.charge_y)
267 dz_h2.Fill(qX, corr_flash_z - e.charge_z)
268 dz_prof.Fill(qX, corr_flash_z - e.charge_z)
272 for ib
in list(range(0, profile_bins)):
274 dy_h1.SetBinContent(ibp, dy_prof.GetBinContent(ibp))
275 dy_h1.SetBinError(ibp, dy_prof.GetBinError(ibp))
276 dy_means[int(ib)] = dy_prof.GetBinContent(ibp)
277 dy_spreads[int(ib)] = dy_prof.GetBinError(ibp)
278 dz_h1.SetBinContent(ibp, dz_prof.GetBinContent(ibp))
279 dz_h1.SetBinError(ibp, dz_prof.GetBinError(ibp))
280 dz_means[int(ib)] = dz_prof.GetBinContent(ibp)
281 dz_spreads[int(ib)] = dz_prof.GetBinError(ibp)
282 rr_h1.SetBinContent(ibp, rr_prof.GetBinContent(ibp))
283 rr_h1.SetBinError(ibp, rr_prof.GetBinError(ibp))
284 rr_means[int(ib)] = rr_prof.GetBinContent(ibp)
285 rr_spreads[int(ib)] = rr_prof.GetBinError(ibp)
286 ratio_h1.SetBinContent(ibp, ratio_prof.GetBinContent(ibp))
287 ratio_h1.SetBinError(ibp, ratio_prof.GetBinError(ibp))
288 ratio_means[int(ib)] = ratio_prof.GetBinContent(ibp)
289 ratio_spreads[int(ib)] = ratio_prof.GetBinError(ibp)
291 for e
in nuslice_tree:
292 if e.slices > e.true_nus:
continue
296 isl = int(qX/xbin_width)
298 if dy_spreads[isl] <= 1.e-8:
299 print(
"Warning zero spread.\n",
300 f
"qX: {qX}. isl: {isl}. dy_spreads[isl]: {dy_spreads[isl]} ")
301 dy_spreads[isl] = dy_spreads[isl+1]
302 if dz_spreads[isl] <= 1.e-8:
303 print(
"Warning zero spread.\n",
304 f
"qX: {qX}. isl: {isl}. dz_spreads[isl]: {dz_spreads[isl]} ")
305 dz_spreads[isl] = dz_spreads[isl+1]
306 if rr_spreads[isl] <= 1.e-8:
307 print(
"Warning zero spread.\n",
308 f
"qX: {qX}. isl: {isl}. rr_spreads[isl]: {rr_spreads[isl]} ")
309 rr_spreads[isl] = rr_spreads[isl+1]
310 if ratio_spreads[isl] <= 1.e-8:
311 print(
"Warning zero spread.\n",
312 f
"qX: {qX}. isl: {isl}. ratio_spreads[isl]: {ratio_spreads[isl]} ")
313 ratio_spreads[isl] = ratio_spreads[isl+1]
316 corr_flash_y = e.flash_yb -
y_bias(e.y_skew, hypo_x, y_bias_slope)
317 corr_flash_z = e.flash_zb -
z_bias(e.z_skew, hypo_x, z_bias_slope)
319 score +=
abs(
abs(corr_flash_y-e.charge_y) - dy_means[isl])/dy_spreads[isl]
320 score +=
abs(
abs(corr_flash_z-e.charge_z) - dz_means[isl])/dz_spreads[isl]
321 score +=
abs(e.flash_rr-rr_means[isl])/rr_spreads[isl]
322 if (detector ==
"sbnd" and pset.UseUncoatedPMT)
or \
323 (detector ==
"icarus" and pset.UseOppVolMetric) :
324 score +=
abs(e.flash_ratio-ratio_means[isl])/ratio_spreads[isl]
326 oldunfolded_score_scatter.Fill(qXGl, e.score)
327 unfolded_score_scatter.Fill(qXGl, score)
328 match_score_scatter.Fill(qX, score)
329 match_score_h1.Fill(score)
331 metrics_filename =
'fm_metrics_' + detector +
'.root'
332 hfile = gROOT.FindObject(metrics_filename)
335 hfile = TFile(metrics_filename,
'RECREATE',
336 'Simple flash matching metrics for ' + detector.upper())
349 fitname_suffix = [
"_h",
"_m",
"_l"]
351 for e,suf
in zip(errors, fitname_suffix):
352 yvals = [a + e*b
for a, b
in zip(rr_means, rr_spreads)]
353 graph = TGraph(x_bins,
354 array(
'f', xvals), array(
'f', yvals))
355 name =
"rr_fit" + suf
356 print(
"Fitting: ", name)
357 f = TF1(name, pset.rr_TF1_fit)
360 rr_fit_funcs.append(f)
365 for e,suf
in zip(errors, fitname_suffix):
366 yvals = [a + e*b
for a, b
in zip(ratio_means, ratio_spreads)]
367 graph = TGraph(x_bins,
368 array(
'f', xvals), array(
'f', yvals))
369 name =
"ratio_fit" + suf
370 print(
"Fitting: ", name)
371 f = TF1(name, pset.ratio_TF1_fit)
374 ratio_fit_funcs.append(f)
375 match_score_scatter.Write()
376 oldunfolded_score_scatter.Write()
377 unfolded_score_scatter.Write()
378 match_score_h1.Write()
381 canv = TCanvas(
"canv")
384 crosses = TGraphErrors(x_bins,
385 array(
'f', xvals), array(
'f', dy_means),
386 array(
'f', xerrs), array(
'f', dy_spreads))
387 crosses.SetLineColor(ROOT.kAzure+9)
388 crosses.SetLineWidth(3)
389 crosses.Draw(
"Psame")
394 crosses = TGraphErrors(x_bins,
395 array(
'f', xvals), array(
'f', dz_means),
396 array(
'f', xerrs), array(
'f', dz_spreads))
397 crosses.SetLineColor(ROOT.kAzure+9)
398 crosses.SetLineWidth(3)
399 crosses.Draw(
"Psame")
404 crosses = TGraphErrors(x_bins,
405 array(
'f', xvals), array(
'f', rr_means),
406 array(
'f', xerrs), array(
'f', rr_spreads))
407 crosses.SetLineColor(ROOT.kAzure+9)
408 crosses.SetLineWidth(3)
409 crosses.Draw(
"Psame")
410 for f
in rr_fit_funcs:
411 f.SetLineColor(ROOT.kOrange+7)
412 f.DrawF1(0., drift_distance,
"lsame")
417 crosses = TGraphErrors(x_bins,
418 array(
'f', xvals), array(
'f', ratio_means),
419 array(
'f', xerrs), array(
'f', ratio_spreads))
420 crosses.SetLineColor(ROOT.kAzure+9)
421 crosses.SetLineWidth(3)
422 crosses.Draw(
"Psame")
423 for f
in ratio_fit_funcs:
424 f.SetLineColor(ROOT.kOrange+7)
425 f.DrawF1(0., drift_distance,
"lsame")
426 canv.Print(
"ratio.pdf")
429 oldunfolded_score_scatter.Draw()
430 canv.Print(
"oldunfolded_score_scatter.pdf")
433 unfolded_score_scatter.Draw()
434 canv.Print(
"unfolded_score_scatter.pdf")
437 match_score_scatter.Draw()
438 canv.Print(
"match_score_scatter.pdf")
441 match_score_h1.Draw()
442 canv.Print(
"match_score.pdf")