make_reco_plots.py
Go to the documentation of this file.
1 import os, sys
2 from array import array
3 from ROOT import *
4 from math import *
5 from optparse import OptionParser
6 
7 import root_style
8 
9 gROOT.SetBatch( kTRUE )
10 root_style.SetRootEnv()
11 
12 def SetHistStyle( hist, fill_color, isData ) :
13 
14  root_style.ApplyAxisStyleAll( hist )
15 
16  hist.SetLineColor( kBlack )
17  hist.SetMarkerColor( kBlack )
18  if isData :
19  hist.SetFillColor( 0 )
20  hist.SetFillStyle( 0 )
21  else :
22  hist.SetFillColor( fill_color )
23  hist.SetFillStyle( 1001 )
24 
25 def GetFillColor( mc_comp ) :
26 
27  fill_color = 0
28  if mc_comp == "Data" : fill_color = 0
29  elif "TotMC" in mc_comp : fill_color = kRed
30  elif "QEL" in mc_comp : fill_color = kAzure+7
31  elif "MEC" in mc_comp : fill_color = kOrange-1
32  elif "RES" in mc_comp : fill_color = kGreen+2
33  elif "DIS" in mc_comp : fill_color = kGray
34  elif "Other" in mc_comp : fill_color = kBlack
35 
36  if opts.mec_tune_shift :
37  if "_TotMC" in mc_comp : fill_color = kRed
38  elif "_QEL" in mc_comp : fill_color = kAzure+7
39  elif "_MEC" in mc_comp : fill_color = kOrange-1
40  elif "_RES" in mc_comp : fill_color = kGreen+2
41  elif "_DIS" in mc_comp : fill_color = kGray
42  elif "_Other" in mc_comp : fill_color = kBlack
43 
44  return fill_color
45 
46 
47 def GetLegendTitle( mc_comp, opts ) :
48 
49  title = "Undefined"
50 
51  if "QEL" in mc_comp : title = "QE"
52  elif "MEC" in mc_comp : title = "MEC"
53  elif "RES" in mc_comp : title = "RES"
54  elif "DIS" in mc_comp : title = "DIS"
55  elif "Other" in mc_comp : title = "Other"
56 
57  if mc_comp == "MEC_MINERvA" : title = "MINERvA 2p2h"
58  elif mc_comp == "MEC" : title = "NOvA 2p2h"
59 
60  if opts.mec_tune_shift :
61  title = "Undefined"
62  if "_QEL" in mc_comp : title = "QEL"
63  elif "_MEC" in mc_comp : title = "MEC"
64  elif "_RES" in mc_comp : title = "RES"
65  elif "_DIS" in mc_comp : title = "DIS"
66  elif "_Other" in mc_comp : title = "Other"
67 
68  '''
69  if opts.no_xsec_weights :
70  title = "Undefined"
71  if "MEC" in mc_comp : title = "MEC"
72  elif "QEL" in mc_comp : title = "QEL"
73  elif "RES" in mc_comp : title = "RES"
74  elif "DIS" in mc_comp : title = "DIS"
75  elif "Other" in mc_comp : title = "Other"
76  '''
77 
78  #if mc_comp == "Data" : title = "NOvA ND Data"
79  if mc_comp == "Data" : title = "ND Data"
80 
81  return title
82 
83 
84 def GetPlotTag( opts ) :
85 
86  tag = ""
87 
88  if opts.mec_tune_shift : tag += "_ShiftMEC" + opts.mec_tune_shift
89 
90  if opts.qel_no_rpa : tag += "_QELNoRPA"
91  if opts.res_no_rpa : tag += "_RESNoRPA"
92 
93  if opts.mec_nominal : tag += "_EmpiricalMEC"
94  elif opts.mec_valencia : tag += "_ValenciaMEC"
95  elif opts.mec_minerva : tag += "_MinervaMEC"
96  elif opts.no_mec : tag += "_NoMEC"
97  elif opts.mec_valencia_empirical : tag += "_ValenciaEmpiricalMEC"
98  if opts.interpolate : tag += "_Interpolate"
99 
100  if opts.default_genie : tag = "_DefaultGENIE"
101 
102  return tag
103 
104 def MakePlots( opts ) :
105 
106  beam = opts.beam
107 
108  filename_data = os.path.join( opts.indir, "hists_2018_%s_data_%s.root" % ( beam, opts.period ) )
109  filename_mc = os.path.join( opts.indir, "hists_2018_%s_mc_%s.root" % ( beam, opts.period ) )
110  if opts.interpolate :
111  filename_mc = os.path.join( opts.indir, "hists_2018_%s_mc_%s_interpolate.root" % ( beam, opts.period ) )
112 
113  #print "filename_mc =", filename_mc
114 
115  if not os.path.exists( filename_data ) : sys.exit( "Data file does not exist:\n%s" % filename_data )
116  if not os.path.exists( filename_mc ) : sys.exit( "MC file does not exist:\n%s" % filename_mc )
117 
118  file_data = TFile( filename_data )
119  file_mc = TFile( filename_mc )
120 
121  print "Hist file data:\n", file_data.GetName()
122  print "Hist file MC:\n", file_mc.GetName()
123 
124  tree_pot_data = file_data.Get( "pot" )
125  pot_data = 0.0
126  for i in range( tree_pot_data.GetEntries() ) :
127  tree_pot_data.GetEntry( i )
128  pot_data += tree_pot_data.pot_data
129 
130  tree_pot_mc = file_mc.Get( "pot" )
131  pot_mc = 0.0
132  for i in range( tree_pot_mc.GetEntries() ) :
133  tree_pot_mc.GetEntry( i )
134  pot_mc += tree_pot_mc.pot_mc
135 
136  mc_pot_scale = pot_data / pot_mc
137 
138  print "pot_data = %.2e" % pot_data
139  print "pot_mc = %.2e" % pot_mc
140  print "mc_pot_scale = %.3f" % mc_pot_scale
141 
142  distributions_temp = []
143  distributions_temp.append( "RecoEhad" )
144  distributions_temp.append( "RecoEhadVis" )
145  distributions_temp.append( "RecoEmu" )
146  distributions_temp.append( "RecoEnu" )
147  distributions_temp.append( "RecoQ2" )
148  distributions_temp.append( "RecoQ3" )
149  distributions_temp.append( "RecoW" )
150  distributions_temp.append( "RecoEhadFrac" )
151  distributions_temp.append( "RecoCosThetaMu" )
152  #distributions_temp.append( "RecoQ2_EhadQuantile1" )
153  #distributions_temp.append( "RecoQ2_EhadQuantile2" )
154  #distributions_temp.append( "RecoQ2_EhadQuantile3" )
155  #distributions_temp.append( "RecoQ2_EhadQuantile4" )
156  distributions_temp.append( "PrimTrkLen" )
157 
158  distributions = []
159 
160  for d in distributions_temp :
161  distributions.append( d )
162  distributions.append( "%s_WSelected" %d )
163  for i in range( 4 ) :
164  distributions.append( "%s_EhadQuantile%d" % ( d, i + 1 ) )
165 
166  if opts.mec_tune_shift :
167  distributions = []
168  distributions.append( "RecoEhadVis" )
169  distributions.append( "RecoQ3" )
170 
171 
172  mc_comps = []
173  mc_comps.append( "MEC" )
174  mc_comps.append( "QEL" )
175  mc_comps.append( "RES" )
176  mc_comps.append( "DIS" )
177  mc_comps.append( "Other" )
178 
179  if opts.mec_nominal :
180  mc_comps[ mc_comps.index( "MEC" ) ] = "MEC_Nominal"
181  elif opts.mec_valencia :
182  mc_comps[ mc_comps.index( "MEC" ) ] = "MEC_Valencia"
183  elif opts.mec_minerva :
184  mc_comps[ mc_comps.index( "MEC" ) ] = "MEC_MINERvA"
185  elif opts.no_mec or opts.mec_valencia_empirical :
186  mc_comps.remove( "MEC" )
187 
188  if opts.qel_no_rpa :
189  mc_comps[ mc_comps.index( "QEL" ) ] = "QEL_NoRPA"
190 
191  if opts.res_no_rpa :
192  mc_comps[ mc_comps.index( "RES" ) ] = "RES_Nominal"
193 
194  if opts.default_genie :
195  mc_comps = []
196  mc_comps.append( "QEL_Nominal" )
197  mc_comps.append( "RES_Nominal" )
198  mc_comps.append( "DIS_Nominal" )
199  mc_comps.append( "Other_Nominal" )
200 
201  elif opts.mec_tune_shift == "Up" :
202  mc_comps = []
203  mc_comps.append( "ShiftsQELike_MEC" )
204  mc_comps.append( "ShiftsQELike_QEL" )
205  mc_comps.append( "ShiftsQELike_RES" )
206  mc_comps.append( "ShiftsQELike_DIS" )
207  mc_comps.append( "ShiftsQELike_Other" )
208 
209  elif opts.mec_tune_shift == "Down" :
210  mc_comps = []
211  mc_comps.append( "ShiftsRESLike_MEC" )
212  mc_comps.append( "ShiftsRESLike_QEL" )
213  mc_comps.append( "ShiftsRESLike_RES" )
214  mc_comps.append( "ShiftsRESLike_DIS" )
215  mc_comps.append( "ShiftsRESLike_Other" )
216 
217  print "mc_comps =", mc_comps
218 
219  canvas_dist = TCanvas()
220  canvas_ratio = TCanvas()
221 
222  #text_size = 0.05
223  text_size = 0.045
224  text_size_wsel = 0.044
225  text_font = 42
226  text_angle = 0.0
227  setNDC = True
228  line_spacing = text_size + 0.01
229 
230  xpos_right = 0.6
231  xpos_left = gStyle.GetPadLeftMargin() + text_size
232  ypos_top = 1.0 - gStyle.GetPadTopMargin() - 0.02
233 
234  isRHC = beam.upper() == "RHC"
235 
236  #text_neutrino = "#nu"
237  #if isRHC : text_neutrino = "#bar{#nu}"
238 
239  #text_neutrino = "#nu_{#mu}"
240  #if isRHC : text_neutrino = "#bar{#nu}_{#mu}"
241  #text_beam = text_neutrino + " Enhanced Beam"
242 
243  text_beam = "Neutrino Beam"
244  if isRHC : text_beam = "Antineutrino Beam"
245  label_beam = root_style.GetLabel( text_beam, xpos_right, ypos_top, text_size, kBlack, text_font, 13, text_angle, setNDC )
246 
247  text_select = "#nu_{#mu} + #bar{#nu}_{#mu} CC Selection"
248  ypos_select = ypos_top - text_size #line_spacing
249  label_select = root_style.GetLabel( text_select, xpos_right, ypos_select, text_size, kBlack, text_font, 13, text_angle, setNDC )
250 
251  label_prelim = root_style.GetLabelPreliminary()
252  label_prelim.SetTextFont( text_font )
253  label_prelim.SetTextSize(2/30.)
254 
255  #pot_power = 20
256  #text_pot = "%.2f #times 10^{%d} P.O.T." % ( pot_data / pow( 10.0, pot_power ), pot_power )
257  #label_pot = root_style.GetLabel( text_pot , xpos_pot , ypos_pot , text_size, kBlack, text_font, 13, text_angle, setNDC )
258 
259  text_list_weights = []
260  if opts.mec_nominal :
261  text_list_weights.append( "Dytman" )
262  text_list_weights.append( "Empirical MEC" )
263  elif opts.mec_valencia :
264  text_list_weights.append( "Valencia MEC" )
265  elif opts.mec_minerva :
266  text_list_weights.append( "MINERvA MEC" )
267  if opts.qel_no_rpa :
268  text_list_weights.append( "No QE RPA" )
269  elif opts.res_no_rpa :
270  text_list_weights.append( "No RES Supp." )
271  elif opts.default_genie :
272  text_list_weights.append( "Default GENIE" )
273 
274  event_scale_power = 4
275  event_scale = 1.0 / pow( 10.0, event_scale_power )
276 
277  for d in distributions :
278 
279  hist_name_data = "%s_Data" % d
280  if not file_data.FindKey( hist_name_data ) :
281  print "Histogram \"%s\" not found. Skipping..." % hist_name_data
282  continue
283  hist_data = file_data.Get( hist_name_data )
284  hist_data.Scale( event_scale )
285  hist_data.GetYaxis().SetTitle( "10^{%d} Events" % event_scale_power )
286  SetHistStyle( hist_data, GetFillColor( "Data" ), True )
287 
288  hists_mc = {}
289  hist_tot_mc_sum = None
290 
291  legend = root_style.GetLegend()
292  legend.SetTextSize( text_size )
293  legend.SetTextFont( text_font )
294  #legend.SetHeader( text_select )
295  legend.AddEntry( hist_data, GetLegendTitle( "Data", opts ), "LPE" )
296 
297  legend_r = root_style.GetLegend()
298  legend_r.SetTextSize( text_size )
299  legend_r.SetTextFont( text_font )
300  legend_r.AddEntry( hist_data, " ", "LPE" )
301 
302  legend_l = root_style.GetLegend()
303  legend_l.SetTextSize( text_size )
304  legend_l.SetTextFont( text_font )
305  legend_l.AddEntry( None, GetLegendTitle( "Data", opts ), "" )
306 
307 
308  for c in mc_comps :
309  hist_name = "%s_%s" % ( d, c )
310  #print "hist_name =", hist_name
311  if not file_mc.FindKey( hist_name ) :
312  print "Histogram \"%s\" not found. Skipping..." % hist_name
313  break
314  hists_mc[ c ] = file_mc.Get( hist_name )
315  hists_mc[ c ].Scale( mc_pot_scale * event_scale )
316  SetHistStyle( hists_mc[ c ], GetFillColor( c ), False )
317 
318  if hist_tot_mc_sum == None :
319  hist_tot_mc_sum = hists_mc[ c ].Clone( "TotMC" )
320  SetHistStyle( hist_tot_mc_sum, GetFillColor( "TotMC" ), False )
321  else :
322  hist_tot_mc_sum.Add( hists_mc[ c ] )
323 
324  #print "len( mc_comps ) =", len( mc_comps )
325  #print "len( hists_mc ) =", len( hists_mc )
326  if len( hists_mc ) != len( mc_comps ) : continue
327 
328  if opts.mec_valencia_empirical :
329 
330  hist_mec_empiricial = file_mc.Get( "%s_MEC_Nominal" % d )
331  hist_mec_empiricial.Scale( mc_pot_scale * event_scale )
332  hist_tot_mc_empirical = hist_tot_mc_sum.Clone( "TotMC_ValenciaMEC" )
333  hist_tot_mc_empirical.Add( hist_mec_empiricial )
334  hist_tot_mc_empirical.SetFillColor( 0 )
335  hist_tot_mc_empirical.SetFillStyle( 0 )
336  hist_tot_mc_empirical.SetLineColor( kRed )
337  #hist_tot_mc_empirical.SetLineColor( kBlue )
338  hist_tot_mc_empirical.SetLineWidth( 3 )
339  legend.AddEntry( hist_tot_mc_empirical, "Empirical MEC", "L" )
340 
341  hist_mec_valencia = file_mc.Get( "%s_MEC_Valencia_FixNP" % d )
342  hist_mec_valencia.Scale( mc_pot_scale * event_scale )
343  hist_tot_mc_valencia = hist_tot_mc_sum.Clone( "TotMC_ValenciaMEC" )
344  hist_tot_mc_valencia.Add( hist_mec_valencia )
345  hist_tot_mc_valencia.SetFillColor( 0 )
346  hist_tot_mc_valencia.SetFillStyle( 0 )
347  #hist_tot_mc_valencia.SetLineColor( kRed )
348  hist_tot_mc_valencia.SetLineColor( kBlack )
349  hist_tot_mc_valencia.SetLineWidth( 3 )
350  hist_tot_mc_valencia.SetLineStyle( 2 )
351  legend.AddEntry( hist_tot_mc_valencia, "Valencia MEC", "L" )
352 
353  for c in mc_comps :
354  legend_name = GetLegendTitle( c, opts )
355  legend.AddEntry( hists_mc[ c ], legend_name, "F" )
356  legend_l.AddEntry( None, legend_name, "" )
357  legend_r.AddEntry( hists_mc[ c ], " ", "F" )
358 
359  stackMC = THStack()
360  for c in reversed( mc_comps ) : stackMC.Add( hists_mc[ c ] )
361 
362  maxData = hist_data.GetBinContent( hist_data.GetMaximumBin() )
363  ymax_scale = 1.5
364  if "RecoCosThetaMu" in d :
365  ymax_scale = 1.3
366  elif "RecoQ2_EhadQuantile3" in d :
367  #if isRHC : ymax_scale = 2.5
368  #else : ymax_scale = 2.0
369  ymax_scale = 2.0
370  elif "RecoQ2_EhadQuantile4" in d :
371  if isRHC : ymax_scale = 2.5
372  else : ymax_scale = 3.0
373  elif d == "RecoEhadVis_EhadQuantile4" :
374  if isRHC : ymax_scale = 2.5
375  else : ymax_scale = 3.0
376  elif d == "RecoQ3_EhadQuantile4" :
377  if isRHC : ymax_scale = 3.5
378  else : ymax_scale = 2.5
379  elif d == "RecoQ2_WSelected" :
380  if isRHC : ymax_scale = 1.9
381  else : ymax_scale = 1.7
382  elif d == "RecoEhadVis" :
383  ymax_scale = 1.3
384  elif d == "RecoEhadFrac" :
385  ymax_scale = 1.3
386 
387  ymax = ymax_scale * maxData
388  stackMC.SetMaximum( ymax )
389  stackMC.SetMinimum( 0.0 )
390  hist_data.SetMaximum( ymax )
391  hist_data.SetMinimum( 0.0 )
392 
393  xmin = hist_data.GetXaxis().GetXmin()
394  xmax = hist_data.GetXaxis().GetXmax()
395  xndiv = hist_data.GetXaxis().GetNdivisions()
396  if "RecoW" in d:
397  xmin = 0.5
398  xmax = 2.5
399  xndiv = 504
400  if "RecoEnu" in d:
401  xmin = 0.0
402  xmax = 4.0
403  xndiv = 508
404  elif "PrimTrkLen" in d :
405  xmax = 15.0
406  #xndiv = 407
407  elif "RecoCosThetaMu" in d :
408  xmin = 0.7
409  xndiv = 506
410  elif "RecoQ2" in d :
411  xmax = 1.0
412  elif "RecoEhadVis" in d :
413  if isRHC :
414  xmax = 0.4
415  xndiv = 504
416  else :
417  xmax = 0.6
418  xndiv = 506
419 
420  if "EhadQuantile" in d :
421  text_qntl = "Quantile %d" % int( d[-1] )
422  label_qntl = root_style.GetLabel( "Quantile X", xpos_left, ypos_top, text_size, kBlack, text_font, 13, text_angle, setNDC )
423 
424  # Set X and Y positions of labels and legend
425 
426  xpos_legend = 0.55
427  width_legend = 1.0 - gStyle.GetPadRightMargin() - 0.005 - xpos_legend
428  xpos_other = gStyle.GetPadLeftMargin() + 0.04
429  if "RecoCosThetaMu" in d or d == "RecoQ3_EhadQuantile4" or d == "RecoW_EhadQuantile4" or "RecoEnu" in d :
430  xpos_legend = gStyle.GetPadLeftMargin() + 0.03 # text_size
431  xpos_other = 0.6
432 
433  label_beam .SetX( xpos_legend + 0.01 )
434  label_select.SetX( xpos_legend + 0.01 )
435  legend.SetX1( xpos_legend )
436  legend.SetX2( xpos_legend + width_legend )
437  legend.SetY2( ypos_select - text_size - 0.02 )
438  legend.SetY1( legend.GetY2() - ( legend.GetTextSize() + 0.01 ) * legend.GetNRows() )
439  legend.SetFillColorAlpha( 0, 0 )
440 
441  '''
442  print "legend.GetTextAlign() =", legend.GetTextAlign()
443 
444  legend_r.SetX1( legend.GetX2() - 0.08 )
445  legend_r.SetX2( legend.GetX2() + width_legend - 0.08 )
446  #legend_r.SetY2( legend.GetY2() +0.003 )
447  #legend_r.SetY1( legend.GetY1() +0.003 )
448  legend_r.SetY2( legend.GetY2() )
449  legend_r.SetY1( legend.GetY1() )
450 
451  legend_r.SetFillColor(kWhite)
452  legend_r.SetFillStyle(0)
453  legend_r.SetLineColor(kWhite)
454 
455  legend_l.SetX1( legend.GetX1() )
456  legend_l.SetX2( legend.GetX2() - 0.1)
457  legend_l.SetY2( legend.GetY2() )
458  legend_l.SetY1( legend.GetY1() )
459  legend_l.SetTextAlign(32)
460  legend_l.SetFillStyle(0)
461  '''
462 
463  '''
464  if d == "RecoQ3_EhadQuantile4" :
465  print "Trying to make legend fill transparent for %s" % d
466  legend.SetFillColorAlpha( 0, 0 )
467  '''
468 
469  text_list_others = []
470  if "EhadQuantile" in d :
471  text_list_others.append( "Quantile %d" % int( d[-1] ) )
472 
473  if "WSelected" in d :
474  text_list_others.append( "1.2 < W_{reco} < 1.5 GeV/c^{2}" )
475 
476  text_list_others += text_list_weights
477  if opts.interpolate :
478  text_list_others.append( "Interpolated MEC" )
479 
480  labels_others = []
481  ypos_other = ypos_top
482  for text in text_list_others :
483  '''
484  if "RecoW" in text :
485  labels_others.append( root_style.GetLabel( text, xpos_other, ypos_other, text_size_wsel, kBlack, text_font, 13, text_angle, setNDC ) )
486  elif "RPA" in text :
487  labels_others.append( root_style.GetLabel( text, xpos_other+.03, ypos_other, text_size_wsel, kBlack, text_font, 13, text_angle, setNDC ) )
488  else :
489  labels_others.append( root_style.GetLabel( text, xpos_other, ypos_other, text_size, kBlack, text_font, 13, text_angle, setNDC ) )
490  '''
491  labels_others.append( root_style.GetLabel( text, xpos_other, ypos_other, text_size, kBlack, text_font, 13, text_angle, setNDC ) )
492  ypos_other -= line_spacing
493 
494  if "RecoQ2" in d :
495  hist_data.GetXaxis().SetTitle( "Reco Q^{2} (GeV^{2}/c^{2})" )
496  elif "RecoQ3" in d :
497  hist_data.GetXaxis().SetTitle( "Reco #left|#vec{q}#right| (GeV/c)" )
498  elif "PrimTrkLen" in d :
499  hist_data.GetXaxis().SetTitle( "Muon Track Length (m)" )
500  elif "RecoW" in d :
501  hist_data.GetXaxis().SetTitle( "Reco W (GeV/c^{2})" )
502 
503  canvas_dist.cd()
504  canvas_dist.Clear()
505 
506  hist_data.GetXaxis().SetRangeUser( xmin, xmax )
507  if "PrimTrkLen" not in d : hist_data.GetXaxis().SetNdivisions( xndiv, False )
508 
509  hist_data.Draw( "E" )
510  if opts.mec_valencia_empirical :
511  hist_tot_mc_empirical.Draw( "hist ][ same" )
512  hist_tot_mc_valencia.Draw( "hist ][ same" )
513  stackMC.Draw( "hist same" )
514  stackMC.GetXaxis().SetRangeUser( xmin, xmax )
515  if opts.draw_tot_mc :
516  hist_tot_mc = file_mc.Get( "%s_TotMC" % d )
517  hist_tot_mc.Scale( mc_pot_scale )
518  hist_tot_mc.SetLineColor( kRed )
519  legend.AddEntry( hist_tot_mc, "Total MC", "L" )
520  hist_tot_mc.Draw( "same hist" )
521  hist_data.Draw( "same E" )
522  if "PrimTrkLen" in d or "RecoEnu" in d :
523  #legend_r.Draw()
524  #legend_l.Draw()
525  legend.Draw()
526  else :
527  legend.Draw()
528  #label_prelim.Draw()
529  label_beam.Draw()
530  label_select.Draw()
531  for label in labels_others : label.Draw()
532 
533  hist_ratio = hist_data.Clone( "%s_Ratio" % d )
534  hist_ratio.Divide( hist_tot_mc_sum )
535  hist_ratio.GetYaxis().SetTitle( "Data / MC" )
536  hist_ratio.GetYaxis().SetRangeUser( 0.5, 1.5 )
537  hist_ratio.GetYaxis().SetNdivisions( 504 )
538 
539  canvas_ratio.cd()
540  canvas_ratio.Clear()
541  canvas_ratio.SetGridy(1)
542  hist_ratio.Draw( "E" )
543 
544  period_tag = opts.period[0].upper() + opts.period[1:]
545 
546  plot_tag = GetPlotTag( opts )
547  if opts.add_out_tag != None : plot_tag += "_%s" % opts.add_out_tag
548 
549  plots_dir = os.path.join( opts.outdir, "%s_%s_Reco%s" % ( beam.upper(), period_tag, plot_tag ) )
550  if not os.path.isdir( plots_dir ) : os.makedirs( plots_dir )
551 
552  out_dist_name = "%s_%s_%s" % ( beam.upper(), period_tag, d )
553  out_dist_name += plot_tag
554 
555  #canvas_dist.SaveAs( "plots/%s.pdf" % out_dist_name )
556  #canvas_dist.SaveAs( os.path.join( plots_dir, "%s.pdf" % out_dist_name ) )
557 
558  root_style.SavePlots( canvas_dist, os.path.join( plots_dir, out_dist_name ) )
559 
560 
561  out_name_ratio = "%s_%s_%s_Ratio" % ( beam.upper(), period_tag, d )
562  out_name_ratio += plot_tag
563 
564  #canvas_ratio.SaveAs( "plots/%s.pdf" % out_name_ratio )
565  #canvas_ratio.SaveAs( os.path.join( plots_dir, "%s.pdf" % out_name_ratio ) )
566 
567 
568 if __name__ == "__main__" :
569 
570  ###############################
571  # Parse Options
572  ###############################
573  #default_indir = "./"
574  default_indir = "/pnfs/nova/persistent/users/mislivec/xsec_tuning_paper/"
575  default_outdir = "./plots/reco/"
576  default_period = "full"
577  help_indir = "Input histogram files directory. Default:\n%s" % default_indir
578  help_outdir = "Output histogram files directory. Default:\n%s" % default_outdir
579  parser = OptionParser()
580  parser.add_option( "--beam" , dest = "beam" , type = "str", default = None, help = "Beam (FHC or RHC)" )
581  parser.add_option( "--period" , dest = "period" , type = "str", default = default_period, help = "Period (default = %s)" % default_period )
582  parser.add_option( "--indir" , dest = "indir" , type = "str", default = default_indir , help = help_indir )
583  parser.add_option( "--outdir" , dest = "outdir" , type = "str", default = default_outdir, help = help_outdir )
584  parser.add_option( "--add_out_tag" , dest = "add_out_tag" , type = "str", default = None, help = "Add tag to output plot file names" )
585  parser.add_option( "--plot_title_all" , dest = "plot_title_all" , type = "str", default = None, help = "Title applied to all plots" )
586  parser.add_option( "--interpolate", dest = "interpolate", action = "store_true", default = False, help = "Interpolated MEC weights" )
587 
588  parser.add_option( "--mec_tune_shift" , dest = "mec_tune_shift" , type = "str" , default = None , help = "Non-mec shift" )
589  parser.add_option( "--qel_no_rpa" , dest = "qel_no_rpa" , action = "store_true", default = False, help = "QEL w/o RPA" )
590  parser.add_option( "--res_no_rpa" , dest = "res_no_rpa" , action = "store_true", default = False, help = "RES w/o RPA" )
591  parser.add_option( "--mec_nominal" , dest = "mec_nominal" , action = "store_true", default = False, help = "Nominal MEC" )
592  parser.add_option( "--no_mec" , dest = "no_mec" , action = "store_true", default = False, help = "No MEC" )
593  parser.add_option( "--mec_valencia" , dest = "mec_valencia" , action = "store_true", default = False, help = "Valencia MEC" )
594  parser.add_option( "--mec_minerva" , dest = "mec_minerva" , action = "store_true", default = False, help = "MINERvA MEC" )
595  parser.add_option( "--draw_tot_mc" , dest = "draw_tot_mc" , action = "store_true", default = False, help = "Draw total MC" )
596  parser.add_option( "--default_genie" , dest = "default_genie" , action = "store_true", default = False, help = "Default GENIE" )
597 
598  parser.add_option( "--mec_valencia_empirical", dest = "mec_valencia_empirical", action = "store_true", default = False, help = "Plot both Valencia and Empirical MEC" )
599 
600  if len(sys.argv) < 2 : parser.parse_args( "--help".split() )
601 
602  ( opts, args ) = parser.parse_args()
603 
604  if not opts.beam : sys.exit( "Beam (FHC or RHC) requred" )
605 
606  MakePlots( opts )
void split(double tt, double *fr)
def SetHistStyle(hist, fill_color, isData)
constexpr T pow(T x)
Definition: pow.h:75
gargamelle SetTitle("Gargamelle #nu_{e} CC data")
def GetPlotTag(opts)
def GetLegendTitle(mc_comp, opts)
simulatedPeak Scale(1/simulationNormalisationFactor)
def GetFillColor(mc_comp)