All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Classes | Functions | Variables
generate_simple_weighted_template Namespace Reference

Classes

class  dotDict
 

Functions

def hypo_flashx_from_H2
 
def x_estimate_and_rms
 
def y_bias
 
def z_bias
 
def generator
 
def main
 

Variables

 myargv = sys.argv
 
string detector = "experiment"
 

Function Documentation

def generate_simple_weighted_template.generator (   nuslice_tree,
  rootfile,
  pset 
)

Definition at line 134 of file generate_simple_weighted_template.py.

135 def generator(nuslice_tree, rootfile, pset):
136  drift_distance = pset.DriftDistance
137  x_bins = pset.XBins
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
142 
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":
149  x_gl_low = -215
150  x_gl_up = 215
151  elif detector == "icarus":
152  if pset.Cryostat == 0:
153  x_gl_low = -380
154  x_gl_up = -50
155  if pset.Cryostat == 1:
156  x_gl_low = 380
157  x_gl_up = 50
158 
159  profile_bins = x_bins
160  profile_option = 's' # errors are the standard deviation
161 
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)")
178 
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)")
195 
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)")
212 
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")
229 
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)")
248 
249  # fill rr_h2 and ratio_h2 first
250  for e in nuslice_tree:
251  if e.slices > e.true_nus: continue
252  qX = e.charge_x
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)
257  # use rr_h2 and ratio_h2 to compute hypo_x
258  for e in nuslice_tree:
259  if e.slices > e.true_nus: continue
260  qX = e.charge_x
261  hypo_x = hypo_flashx_from_H2(e.flash_rr, rr_h2, e.flash_ratio, ratio_h2)
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)
264 
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)
269 
270 
271  # fill histograms for match score calculation from profile histograms
272  for ib in list(range(0, profile_bins)):
273  ibp = ib + 1
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)
290 
291  for e in nuslice_tree:
292  if e.slices > e.true_nus: continue
293  qX = e.charge_x
294  qXGl = e.charge_x_gl
295  # calculate match score
296  isl = int(qX/xbin_width)
297  score = 0.
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]
314  # ratio_spreads[isl] = ratio_spreads[isl-1]
315  hypo_x = hypo_flashx_from_H2(e.flash_rr, rr_h2, e.flash_ratio, ratio_h2)
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)
318 
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]
325 
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)
330 
331  metrics_filename = 'fm_metrics_' + detector + '.root'
332  hfile = gROOT.FindObject(metrics_filename)
333  if hfile:
334  hfile.Close()
335  hfile = TFile(metrics_filename, 'RECREATE',
336  'Simple flash matching metrics for ' + detector.upper())
337 
338  dy_h2.Write()
339  dy_prof.Write()
340  dy_h1.Write()
341  dz_h2.Write()
342  dz_prof.Write()
343  dz_h1.Write()
344  rr_h2.Write()
345  rr_prof.Write()
346  rr_h1.Write()
347 
348  errors = [1, 0, -1]
349  fitname_suffix = ["_h", "_m", "_l"]
350  rr_fit_funcs = []
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)
358  graph.Fit(f)
359  f.Write()
360  rr_fit_funcs.append(f)
361  ratio_h2.Write()
362  ratio_prof.Write()
363  ratio_h1.Write()
364  ratio_fit_funcs = []
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)
372  graph.Fit(f)
373  f.Write()
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()
379  hfile.Close()
380 
381  canv = TCanvas("canv")
382 
383  dy_h2.Draw()
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")
390  canv.Print("dy.pdf")
391  canv.Update()
392 
393  dz_h2.Draw()
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")
400  canv.Print("dz.pdf")
401  canv.Update()
402 
403  rr_h2.Draw()
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")
413  canv.Print("rr.pdf")
414  canv.Update()
415 
416  ratio_h2.Draw()
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")
427  canv.Update()
428 
429  oldunfolded_score_scatter.Draw()
430  canv.Print("oldunfolded_score_scatter.pdf")
431  canv.Update()
432 
433  unfolded_score_scatter.Draw()
434  canv.Print("unfolded_score_scatter.pdf")
435  canv.Update()
436 
437  match_score_scatter.Draw()
438  canv.Print("match_score_scatter.pdf")
439  canv.Update()
440 
441  match_score_h1.Draw()
442  canv.Print("match_score.pdf")
443  canv.Update()
444  sleep(20)
445 
# Main program.
do one_file $F done echo for F in find $TOP name CMakeLists txt print
T abs(T value)
auto zip(Iterables &&...iterables)
Range-for loop helper iterating across many collections at the same time.
Definition: zip.h:295
list
Definition: file_to_url.sh:28
def generate_simple_weighted_template.hypo_flashx_from_H2 (   flash_rr,
  rr_h2,
  flash_ratio,
  ratio_h2 
)

Definition at line 83 of file generate_simple_weighted_template.py.

83 
84 def hypo_flashx_from_H2(flash_rr, rr_h2, flash_ratio, ratio_h2):
85  rr_hypoX, rr_hypoXWgt = x_estimate_and_rms(flash_rr, rr_h2);
86  ratio_hypoX, ratio_hypoXWgt = x_estimate_and_rms(flash_ratio, ratio_h2);
87 
88  sum_weights = rr_hypoXWgt + ratio_hypoXWgt
89  hypo_x = (rr_hypoX*rr_hypoXWgt + ratio_hypoX*ratio_hypoXWgt) / sum_weights
90  hypo_x_err = np.sqrt(sum_weights) / sum_weights
91  # return (hypo_x, hypo_x_err, rr_hypoX, ratio_hypoX)
92  return hypo_x
93 
def generate_simple_weighted_template.main ( )

Definition at line 446 of file generate_simple_weighted_template.py.

447 def main():
448 
449  # Parse arguments.
450  parser = argparse.ArgumentParser(prog='generate_simple_weighted_template.py')
451  parser.add_argument('file')
452  # parser.add_argument('--help', '-h',
453  # action='store_true',
454  # help='help flag' )
455  group = parser.add_mutually_exclusive_group(required=True)
456  group.add_argument('--sbnd', action='store_true', help='Generate metrics for SBND')
457  group.add_argument('--icarus', action='store_true', help='Generate metrics for ICARUS')
458  args = parser.parse_args()
459 
460  # if args.help:
461  # print("To run do:\n"/
462  # "generate_simple_weighted_template.py file.root\n"/
463  # "where file.root has a fmatch/nuslicetree")
464  # return(0)
465  if args.sbnd :
466  print("Generate metrics for SBND")
467  elif args.icarus :
468  print("Generate metrics for ICARUS")
469 
470  if not larbatch_posix.exists(args.file):
471  print('Input file %s does not exist.' % args.file)
472  return 1
473 
474  print('\nOpening %s' % args.file)
475  rootfile = TFile.Open(args.file)
476  if not rootfile.IsOpen() or rootfile.IsZombie():
477  print('Failed to open %s' % args.file)
478  return 1
479  global detector
480  if args.sbnd:
481  fcl_params = fhicl.make_pset('flashmatch_sbnd.fcl')
482  pset = dotDict(fcl_params['sbnd_simple_flashmatch'])
483  detector = "sbnd"
484  dir = rootfile.Get(args.file+":/fmatch")
485  nuslice_tree = dir.Get("nuslicetree") # , nuslice_tree)
486  # nuslice_tree.Print()
487  elif args.icarus:
488  fcl_params = fhicl.make_pset('flashmatch_simple_icarus.fcl')
489  # TODO: add option to use cryo 0 and cryo 1
490  pset = dotDict(fcl_params['icarus_simple_flashmatch_0'])
491  detector = "icarus"
492  dir = rootfile.Get(args.file+":/fmatchCryo0")
493  nuslice_tree = dir.Get("nuslicetree") # , nuslice_tree)
494  # nuslice_tree.Print()
495 
496  generator(nuslice_tree, rootfile, pset)
497 
do one_file $F done echo for F in find $TOP name CMakeLists txt print
def generate_simple_weighted_template.x_estimate_and_rms (   metric_value,
  metric_h2 
)

Definition at line 94 of file generate_simple_weighted_template.py.

94 
95 def x_estimate_and_rms(metric_value, metric_h2):
96  kMinEntriesInProjection = 100
97  fXBinWidth = 5.
98  bin = metric_h2.GetYaxis().FindBin(metric_value);
99  bins = metric_h2.GetNbinsY();
100  metric_hypoX = -1.;
101  metric_hypoXWgt = 0.;
102  bin_buff = 0;
103  while 0 < bin-bin_buff or bin+bin_buff <= bins :
104  low_bin = bin-bin_buff if 0 < bin-bin_buff else 0
105  high_bin = bin+bin_buff if bin+bin_buff <= bins else -1
106  metric_px = metric_h2.ProjectionX("metric_px", low_bin, high_bin);
107  if metric_px.GetEntries() > kMinEntriesInProjection :
108  metric_hypoX = metric_px.GetRandom();
109  metric_rmsX = metric_px.GetRMS();
110  if metric_rmsX < fXBinWidth: # something went wrong
111  print(f"metric_h2 projected on metric_value: {metric_value}, bin: {bin}, "
112  f"bin_buff: {bin_buff}; has {metric_px.GetEntries()} entries.")
113  print(f"metric_hypoX: {metric_hypoX}, metric_rmsX: {metric_rmsX}")
114  return (-1., 0.); # no estimate
115  metric_hypoXWgt = 1/(metric_rmsX*metric_rmsX);
116  return (metric_hypoX, metric_hypoXWgt);
117  bin_buff += 1;
118  return (-1., 0.); # no estimate
119 
do one_file $F done echo for F in find $TOP name CMakeLists txt print
def generate_simple_weighted_template.y_bias (   y_skew,
  hypo_x,
  y_bias_slope 
)

Definition at line 120 of file generate_simple_weighted_template.py.

121 def y_bias(y_skew, hypo_x, y_bias_slope):
122  if detector == "sbnd":
123  return y_skew * hypo_x * hypo_x * y_bias_slope
124  elif detector == "icarus":
125  return y_skew * hypo_x * y_bias_slope
126 
def generate_simple_weighted_template.z_bias (   z_skew,
  hypo_x,
  z_bias_slope 
)

Definition at line 127 of file generate_simple_weighted_template.py.

128 def z_bias(z_skew, hypo_x, z_bias_slope):
129  if detector == "sbnd":
130  return z_skew * hypo_x * z_bias_slope;
131  elif detector == "icarus":
132  return z_skew * hypo_x * z_bias_slope;
133 

Variable Documentation

string generate_simple_weighted_template.detector = "experiment"

Definition at line 61 of file generate_simple_weighted_template.py.

generate_simple_weighted_template.myargv = sys.argv

Definition at line 54 of file generate_simple_weighted_template.py.