make_true_q0q3_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 #import canMan
9 
10 gROOT.SetBatch( kTRUE )
11 root_style.SetRootEnv()
12 
13 def RebinHistQ0Q3( hist ) :
14 
15  hist.RebinX( 2 )
16  hist.RebinY( 2 )
17 
18 def SetRangeQ0Q3( hist ) :
19 
20  hist.GetXaxis().SetRangeUser( 0.0, 1.0 )
21  hist.GetYaxis().SetRangeUser( 0.0, 0.8 )
22  #hist.GetYaxis().SetLimits( 1.0e-3, 0.8 / hist.GetYaxis().GetXmax() )
23  #print hist.GetYaxis().GetXmax()
24 
25 def GetMaximum( hist ) :
26 
27  return hist.GetBinContent( hist.GetMaximumBin() )
28 
30 
31  # A colour palette that goes blue->white->red, useful for correlation matrices
32 
33  n_color_contours = 999
34  #n_color_contours = 25
35  NRGBs = 3;
36  initialized = False
37  colors = array( "i" )
38 
39  stops = array( "d", [ 0.00, 0.20, 1.00 ] )
40  red = array( "d", [ 0.00, 1.00, 1.00 ] )
41  green = array( "d", [ 0.00, 1.00, 0.00 ] )
42  blue = array( "d", [ 1.00, 1.00, 0.00 ] )
43  colmin = TColor.CreateGradientColorTable( NRGBs, stops, red, green, blue, n_color_contours )
44  for i in range( n_color_contours ) : colors.insert( i, colmin + i )
45 
46  gStyle.SetNumberContours( n_color_contours )
47  gStyle.SetPalette( n_color_contours, colors )
48 
49 def SetTitleQ0( hist ) :
50 
51  hist.GetXaxis().SetTitle( "True #left|#vec{q}#right| (GeV/c)" )
52 
53 def SetDecimalsAllAxes( hist ) :
54 
55  hist.GetXaxis().SetDecimals( kTRUE )
56  hist.GetYaxis().SetDecimals( kTRUE )
57  hist.GetZaxis().SetDecimals( kTRUE )
58 
59 def SetHistLineAtt( hist, width, color ) :
60 
61  hist.SetLineWidth( width )
62  hist.SetLineColor( color )
63 
64 if __name__ == "__main__":
65 
66  ###############################
67  # Parse Options
68  ###############################
69 
70  #default_indir = "./"
71  default_indir = "/pnfs/nova/persistent/users/mislivec/xsec_tuning_paper/"
72  default_outdir = './plots/true_q0q3'
73 
74  #default_indir = "/pnfs/nova/persistent/users/mislivec/xs_tuning_hists/rpa_fix"
75  #default_outdir = "/nova/app/users/mislivec/scripts_xs_tuning_blessed_plots/plots_true_q0q3_rpa_fix"
76 
77  parser = OptionParser()
78  parser.add_option( '--beam' , dest = 'beam' , default = None , help = 'Beam (FHC or RHC)' )
79  parser.add_option( '--sim', dest = 'sim', action = "store_true", default = False, help = "Add \"NOvA Simulation\" label" )
80  #parser.add_option( '--outdir', dest = 'outdir', default = default_outdir, help = 'Output plots directory:\n%s' % default_outdir )
81  #parser.add_option( '--indir' , dest = 'indir' , default = default_indir , help = 'Input plots directory:\n%s' % default_indir )
82 
83  if len(sys.argv) < 2 : parser.parse_args( '--help'.split() )
84 
85  ( opts, args ) = parser.parse_args()
86 
87  if not opts.beam : sys.exit( 'Beam (FHC or RHC) requred' )
88 
89  #hists_file_name_data = os.path.join( opts.indir, "hists_2018_%s_data_full.root" % opts.beam )
90  #hists_file_name_mc = os.path.join( opts.indir, "hists_2018_%s_mc_full.root" % opts.beam )
91 
92  hists_file_name_data = os.path.join( default_indir, "hists_2018_%s_data_full.root" % opts.beam )
93  hists_file_name_mc = os.path.join( default_indir, "hists_2018_%s_mc_full.root" % opts.beam )
94 
95  if not os.path.exists( hists_file_name_data ) : sys.exit( 'File does not exist:\n%s' % hists_file_name_data )
96  if not os.path.exists( hists_file_name_mc ) : sys.exit( 'File does not exist:\n%s' % hists_file_name_mc )
97  print 'Input MC histogram file:\n', hists_file_name_data
98  print 'Input MC histogram file:\n', hists_file_name_mc
99 
100  hists_file_data = TFile( hists_file_name_data )
101  hists_file_mc = TFile( hists_file_name_mc )
102 
103  tree_pot_data = hists_file_data.Get( 'pot' )
104  pot_data = 0.0
105  for i in range( tree_pot_data.GetEntries() ) :
106  tree_pot_data.GetEntry( i )
107  pot_data += tree_pot_data.pot_data
108 
109  tree_pot_mc = hists_file_mc.Get( 'pot' )
110  pot_mc = 0.0
111  for i in range( tree_pot_mc.GetEntries() ) :
112  tree_pot_mc.GetEntry( i )
113  pot_mc += tree_pot_mc.pot_mc
114 
115  mc_pot_scale = pot_data / pot_mc
116 
117  print 'pot_data = %.2e' % pot_data
118  print 'pot_mc = %.2e' % pot_mc
119  print 'mc_pot_scale = %.3f' % mc_pot_scale
120 
121  event_scale_power_1d = 4
122  event_scale_1d = 1.0 / pow( 10.0, event_scale_power_1d )
123  text_events_1d = '10^{%d} Events' % event_scale_power_1d
124  scale_hist_1d = mc_pot_scale * event_scale_1d
125 
126  for k in hists_file_mc.GetListOfKeys() :
127  name = k.GetName()
128  #print name
129  #if 'TrueQ0_vs_TrueQ3' in name : print name
130  #if 'TrueQ0_vs_TrueQ3_MECShapeUp_MEC_' in name : print name
131  #if "TrueQ0_vs_TrueQ3_MECShape" in name : print name
132 
133  isRHC = opts.beam.upper() == 'RHC'
134 
135  tag_neutrino_mec_weights = 'Nu'
136  if isRHC : tag_neutrino_mec_weights = 'Nubar'
137 
138  weights_file_name = os.path.join( os.environ[ 'SRT_PUBLIC_CONTEXT' ], "MCReweight/data/xs/rw_empiricalMEC2018_RPAfix_%s.root" % tag_neutrino_mec_weights.lower() )
139  if not os.path.exists( weights_file_name ) : sys.exit( 'File does not exist:\n%s' % weights_file_name )
140 
141  print 'Input MEC weights file:'
142  print weights_file_name
143 
144  #sys.exit()
145 
146  weights_file = TFile( weights_file_name )
147  for k in weights_file.GetListOfKeys() : print k.GetName()
148 
149  #plots_dir = opts.outdir
150  plots_dir = default_outdir
151  if not os.path.isdir( plots_dir ) : os.makedirs( plots_dir )
152 
153  tag_neutrino_mec_hists = 'Numu'
154  if isRHC : tag_neutrino_mec_hists = 'Numubar'
155 
156  hist_nominal = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MEC_Nominal_%s' % tag_neutrino_mec_hists )
157  hist_tuned = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MEC_%s' % tag_neutrino_mec_hists )
158  hist_valencia = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MEC_Valencia_%s' % tag_neutrino_mec_hists )
159  hist_minerva = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MEC_MINERvA_%s' % tag_neutrino_mec_hists )
160  hist_shift_up = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MECShapeUp_MEC_%s' % tag_neutrino_mec_hists )
161  hist_shift_down = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MECShapeDown_MEC_%s' % tag_neutrino_mec_hists )
162  hist_weights = weights_file .Get( '%s_mec_weights_smoothed' % tag_neutrino_mec_hists.lower() )
163 
164  '''
165  hist_all_mec = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MEC' )
166  hist_all_mec_shift_up = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MECShapeUp_MEC' )
167  hist_all_mec_shift_down = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MECShapeDown_MEC' )
168  hist_non_mec = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_TotMC' )
169  hist_non_mec_shift_up = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MECShapeUp_TotMC' )
170  hist_non_mec_shift_down = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MECShapeDown_TotMC' )
171  hist_non_mec .Add( hist_all_mec, -1.0 )
172  hist_non_mec_shift_up .Add( hist_all_mec_shift_up, -1.0 )
173  hist_non_mec_shift_down.Add( hist_all_mec_shift_down, -1.0 )
174  '''
175 
176  #hist_all_mec = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MEC' )
177  #hist_all_mec_shift_up = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_ShiftsQELike_MEC' )
178  #hist_all_mec_shift_down = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_ShiftsRESLike_MEC' )
179  hist_non_mec = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_TotMC' )
180  hist_non_mec_shift_up = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_ShiftsQELike_TotMC' )
181  hist_non_mec_shift_down = hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_ShiftsRESLike_TotMC' )
182  hist_non_mec .Add( hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_MEC' ), -1.0 )
183  hist_non_mec_shift_up .Add( hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_ShiftsQELike_MEC' ), -1.0 )
184  hist_non_mec_shift_down.Add( hists_file_mc.Get( 'TrueQ0_vs_TrueQ3_ShiftsRESLike_MEC' ), -1.0 )
185 
186 
187  root_style.ApplyAxisStyleAll( hist_nominal )
188  root_style.ApplyAxisStyleAll( hist_tuned )
189  root_style.ApplyAxisStyleAll( hist_valencia )
190  root_style.ApplyAxisStyleAll( hist_minerva )
191  root_style.ApplyAxisStyleAll( hist_shift_up )
192  root_style.ApplyAxisStyleAll( hist_shift_down )
193  root_style.ApplyAxisStyleAll( hist_weights )
194  root_style.ApplyAxisStyleAll( hist_non_mec )
195  root_style.ApplyAxisStyleAll( hist_non_mec_shift_up )
196  root_style.ApplyAxisStyleAll( hist_non_mec_shift_down )
197 
198  '''
199  hist_nominal .Scale( mc_pot_scale * event_scale_1d )
200  hist_tuned .Scale( mc_pot_scale * event_scale_1d )
201  hist_valencia .Scale( mc_pot_scale * event_scale_1d )
202  hist_minerva .Scale( mc_pot_scale * event_scale_1d )
203  hist_shift_up .Scale( mc_pot_scale * event_scale_1d )
204  hist_shift_down .Scale( mc_pot_scale * event_scale_1d )
205  hist_non_mec .Scale( mc_pot_scale * event_scale_1d )
206  hist_non_mec_shift_up .Scale( mc_pot_scale * event_scale_1d )
207  hist_non_mec_shift_down.Scale( mc_pot_scale * event_scale_1d )
208  '''
209 
210  RebinHistQ0Q3( hist_nominal )
211  RebinHistQ0Q3( hist_tuned )
212  RebinHistQ0Q3( hist_valencia )
213  RebinHistQ0Q3( hist_minerva )
214  RebinHistQ0Q3( hist_shift_up )
215  RebinHistQ0Q3( hist_shift_down )
216  RebinHistQ0Q3( hist_non_mec )
217  RebinHistQ0Q3( hist_non_mec_shift_up )
218  RebinHistQ0Q3( hist_non_mec_shift_down )
219 
220  SetTitleQ0( hist_nominal )
221  SetTitleQ0( hist_tuned )
222  SetTitleQ0( hist_valencia )
223  SetTitleQ0( hist_minerva )
224  SetTitleQ0( hist_shift_up )
225  SetTitleQ0( hist_shift_down )
226  SetTitleQ0( hist_non_mec )
227  SetTitleQ0( hist_non_mec_shift_up )
228  SetTitleQ0( hist_non_mec_shift_down )
229  SetTitleQ0( hist_weights )
230 
231  SetDecimalsAllAxes( hist_nominal )
232  SetDecimalsAllAxes( hist_tuned )
233  SetDecimalsAllAxes( hist_weights )
234  #SetDecimalsAllAxes( )
235 
236  margin_left = gStyle.GetPadLeftMargin()
237  margin_right = gStyle.GetPadRightMargin()
238  margin_top = gStyle.GetPadTopMargin()
239  text_size = 0.05
240  text_font = 42
241  text_angle = 0.0
242  setNDC = True
243 
244  xpos_label_mec = 0.37
245  xpos_label_no_mec = 0.37
246  xpos_legend_mec = xpos_label_mec - 0.02
247  xpos_legend_no_mec = xpos_label_no_mec - 0.02
248  ypos_label_top = 1.0 - gStyle.GetPadTopMargin() - 0.02
249  ypos_label_mec = ypos_label_top - text_size
250  ypos_label_no_mec = ypos_label_top - text_size
251  if not isRHC :
252  ypos_label_mec -= 0.005
253  #ypos_label_no_mec -= 0.005
254  ypos_legend_mec = ypos_label_mec - text_size - 0.01
255  ypos_legend_no_mec = ypos_label_no_mec - text_size - 0.02
256 
257  text_mec_neutrino = '#nu_{#mu}'
258  if isRHC : text_mec_neutrino = '#bar{#nu}_{#mu}'
259 
260  text_beam = "Neutrino Beam"
261  if isRHC : text_beam = "Antineutrino Beam"
262 
263  text_nova = "NOvA" #"NO#nuA"
264  text_minerva = "MINERvA" #"MINER#nuA"
265  #text_mec = "%s CC MEC" % text_mec_neutrino
266  text_numu_cc = "%s CC" % text_mec_neutrino
267  text_no_2p2h = "#nu_{#mu} + #bar{#nu}_{#mu} CC No 2p2h"
268  text_empirical_mec = "Empirical MEC"
269  text_valencia_mec = "Valencia MEC"
270  text_nova_tune = "2018 %s #nu + #bar{#nu} Tune" % text_nova
271  text_nova_2p2h = "%s 2p2h" % text_nova
272  text_minerva_2p2h = "%s 2p2h" % text_minerva
273  text_simulation = "%s Simulation" % text_nova
274 
275  label_beam = root_style.GetLabel( text_beam , xpos_label_mec , ypos_label_top , text_size, kBlack, text_font, 13, text_angle, setNDC )
276  #label_mec = root_style.GetLabel( text_mec , xpos_label_mec , ypos_label_mec , text_size, kBlack, text_font, 13, text_angle, setNDC )
277  label_numu_cc = root_style.GetLabel( text_numu_cc, xpos_label_mec , ypos_label_mec , text_size, kBlack, text_font, 13, text_angle, setNDC )
278  label_no_mec = root_style.GetLabel( text_no_2p2h, xpos_label_no_mec, ypos_label_no_mec, text_size, kBlack, text_font, 13, text_angle, setNDC )
279  label_simulation = root_style.GetLabel( text_simulation, 1.0 - margin_right, 1.0 - margin_top * 0.7, 2/30., kBlue , text_font, 31, text_angle, setNDC )
280 
281  canvas = TCanvas()
282 
283  for var in ( 'TrueQ0', 'TrueQ3' ) :
284 
285  xmax = 1.0
286  if var == 'TrueQ0' :
287  hist_1D_nominal = hist_nominal .ProjectionY( '%s_Empirical' % var )
288  hist_1D_tuned = hist_tuned .ProjectionY( '%s_2018CV' % var )
289  hist_1D_valencia = hist_valencia .ProjectionY( '%s_Valencia' % var )
290  hist_1D_minerva = hist_minerva .ProjectionY( '%s_MINERvA' % var )
291  hist_1D_shift_up = hist_shift_up .ProjectionY( '%s_ShiftUp' % var )
292  hist_1D_shift_down = hist_shift_down.ProjectionY( '%s_ShiftDown' % var )
293 
294  hist_1D_non_mec = hist_non_mec .ProjectionY( '%s_NonMEC' % var )
295  hist_1D_non_mec_shift_up = hist_non_mec_shift_up .ProjectionY( '%s_NonMEC_ShiftUp' % var )
296  hist_1D_non_mec_shift_down = hist_non_mec_shift_down.ProjectionY( '%s_NonMEC_ShiftDown' % var )
297 
298  xmax = 1.0
299  else :
300  hist_1D_nominal = hist_nominal .ProjectionX( '%s_Empirical' % var )
301  hist_1D_tuned = hist_tuned .ProjectionX( '%s_2018CV' % var )
302  hist_1D_valencia = hist_valencia .ProjectionX( '%s_Valencia' % var )
303  hist_1D_minerva = hist_minerva .ProjectionX( '%s_MINERvA' % var )
304  hist_1D_shift_up = hist_shift_up .ProjectionX( '%s_ShiftUp' % var )
305  hist_1D_shift_down = hist_shift_down.ProjectionX( '%s_ShiftDown' % var )
306 
307  hist_1D_non_mec = hist_non_mec .ProjectionX( '%s_NonMEC' % var )
308  hist_1D_non_mec_shift_up = hist_non_mec_shift_up .ProjectionX( '%s_NonMEC_ShiftUp' % var )
309  hist_1D_non_mec_shift_down = hist_non_mec_shift_down.ProjectionX( '%s_NonMEC_ShiftDown' % var )
310 
311  xmax = 1.2
312 
313  root_style.ApplyAxisStyleAll( hist_1D_nominal )
314  root_style.ApplyAxisStyleAll( hist_1D_tuned )
315  root_style.ApplyAxisStyleAll( hist_1D_valencia )
316  root_style.ApplyAxisStyleAll( hist_1D_minerva )
317  root_style.ApplyAxisStyleAll( hist_1D_shift_up )
318  root_style.ApplyAxisStyleAll( hist_1D_shift_down )
319  root_style.ApplyAxisStyleAll( hist_1D_non_mec )
320  root_style.ApplyAxisStyleAll( hist_1D_non_mec_shift_up )
321  root_style.ApplyAxisStyleAll( hist_1D_non_mec_shift_down )
322 
323  hist_1D_nominal .GetXaxis().SetRangeUser( 0, xmax )
324  hist_1D_tuned .GetXaxis().SetRangeUser( 0, xmax )
325  hist_1D_valencia .GetXaxis().SetRangeUser( 0, xmax )
326  hist_1D_minerva .GetXaxis().SetRangeUser( 0, xmax )
327  hist_1D_shift_up .GetXaxis().SetRangeUser( 0, xmax )
328  hist_1D_shift_down.GetXaxis().SetRangeUser( 0, xmax )
329  #hist_1D_non_mec .GetXaxis().SetRangeUser( 0, xmax )
330  #hist_1D_non_mec_shift_up .GetXaxis().SetRangeUser( 0, xmax )
331  #hist_1D_non_mec_shift_down.GetXaxis().SetRangeUser( 0, xmax )
332 
333  hist_1D_nominal .Scale( scale_hist_1d )
334  hist_1D_tuned .Scale( scale_hist_1d )
335  hist_1D_valencia .Scale( scale_hist_1d )
336  hist_1D_minerva .Scale( scale_hist_1d )
337  hist_1D_shift_up .Scale( scale_hist_1d )
338  hist_1D_shift_down .Scale( scale_hist_1d )
339  hist_1D_non_mec .Scale( scale_hist_1d )
340  hist_1D_non_mec_shift_up .Scale( scale_hist_1d )
341  hist_1D_non_mec_shift_down.Scale( scale_hist_1d )
342 
343  hist_1D_nominal .GetYaxis().SetTitle( text_events_1d )
344  hist_1D_tuned .GetYaxis().SetTitle( text_events_1d )
345  hist_1D_valencia .GetYaxis().SetTitle( text_events_1d )
346  hist_1D_minerva .GetYaxis().SetTitle( text_events_1d )
347  hist_1D_shift_up .GetYaxis().SetTitle( text_events_1d )
348  hist_1D_shift_down .GetYaxis().SetTitle( text_events_1d )
349  hist_1D_non_mec .GetYaxis().SetTitle( text_events_1d )
350  hist_1D_non_mec_shift_up .GetYaxis().SetTitle( text_events_1d )
351  hist_1D_non_mec_shift_down.GetYaxis().SetTitle( text_events_1d )
352 
353  hist_line_width = 3
354  SetHistLineAtt( hist_1D_nominal , hist_line_width, kGreen+2 )
355  SetHistLineAtt( hist_1D_tuned , hist_line_width, kBlack )
356  SetHistLineAtt( hist_1D_valencia , hist_line_width, kRed )
357  SetHistLineAtt( hist_1D_minerva , hist_line_width, kBlue )
358  SetHistLineAtt( hist_1D_shift_up , hist_line_width, kBlue )
359  SetHistLineAtt( hist_1D_shift_down, hist_line_width, kRed )
360  SetHistLineAtt( hist_1D_non_mec , hist_line_width, kBlack )
361  SetHistLineAtt( hist_1D_non_mec_shift_up , hist_line_width, kBlue )
362  SetHistLineAtt( hist_1D_non_mec_shift_down, hist_line_width, kRed )
363 
364  max_nominal = GetMaximum( hist_1D_nominal )
365  max_tuned = GetMaximum( hist_1D_tuned )
366  max_valencia = GetMaximum( hist_1D_valencia )
367  max_minerva = GetMaximum( hist_1D_minerva )
368  max_shift_up = GetMaximum( hist_1D_shift_up )
369  max_shift_down = GetMaximum( hist_1D_shift_down )
370  max_non_mec = GetMaximum( hist_1D_non_mec )
371  max_non_mec_shift_up = GetMaximum( hist_1D_non_mec_shift_up )
372  max_non_mec_shift_down = GetMaximum( hist_1D_non_mec_shift_down )
373 
374  legend1 = root_style.GetLegend()
375  legend1.SetTextSize( text_size )
376  legend1.SetTextFont( text_font )
377  #legend1.SetHeader( '%s CC MEC' % text_mec_neutrino )
378  legend1.AddEntry( hist_1D_tuned , text_nova_2p2h , 'L' )
379  legend1.AddEntry( hist_1D_nominal, text_empirical_mec, 'L' )
380  legend1.SetX1( xpos_legend_mec )
381  legend1.SetY2( ypos_legend_mec )
382  root_style.SetLegendPosY1( legend1 )
383 
384  canvas.Clear()
385 
386  hist_1D_tuned.SetMaximum( 1.8 * max( max_nominal, max_tuned ) )
387  hist_1D_tuned.Draw( 'hist' )
388  hist_1D_nominal.Draw( 'hist same' )
389  hist_1D_tuned.Draw( 'hist same' )
390  legend1.Draw()
391  if opts.sim : label_simulation.Draw()
392  label_beam.Draw()
393  label_numu_cc.Draw()
394  root_style.SavePlots( canvas, os.path.join( plots_dir, '%s_%s_MEC_2018CV_vs_Empirical' % ( tag_neutrino_mec_hists, var ) ) )
395 
396  text_variable = 'blah'
397  if var == 'TrueQ0' :
398  text_variable = 'true energy transfer, q_0,'
399  else :
400  text_variable = 'true three-momentum transfer magnitude, |q|,'
401 
402  text_caption = 'The %s near detector %s distribution for true %s CC MEC events that pass the 2018 numu selection. ' % ( opts.beam.upper(), text_variable, tag_neutrino_mec_hists.lower() )
403 
404  caption = open( os.path.join( plots_dir, '%s_%s_MEC_2018CV_vs_Empirical.txt' % ( tag_neutrino_mec_hists, var ) ), 'w' )
405  caption.write( text_caption + 'Shown are the 2018 CV MEC and GENIE Empirical MEC.' )
406  caption.close()
407 
408  legend2 = root_style.GetLegend()
409  legend2.SetTextSize( text_size )
410  legend2.SetTextFont( text_font )
411  #legend2.SetHeader( '%s CC MEC' % text_mec_neutrino )
412  legend2.AddEntry( hist_1D_tuned , text_nova_2p2h , 'L' )
413  legend2.AddEntry( hist_1D_valencia, text_valencia_mec, 'L' )
414  legend2.SetX1( xpos_legend_mec )
415  legend2.SetY2( ypos_legend_mec )
416  root_style.SetLegendPosY1( legend2 )
417 
418  canvas.Clear()
419 
420  hist_1D_tuned.SetMaximum( 1.8 * max( max_valencia, max_tuned ) )
421  hist_1D_tuned.Draw( 'hist' )
422  hist_1D_valencia.Draw( 'hist same' )
423  hist_1D_tuned.Draw( 'hist same' )
424  legend2.Draw()
425  if opts.sim : label_simulation.Draw()
426  label_beam.Draw()
427  label_numu_cc.Draw()
428  root_style.SavePlots( canvas, os.path.join( plots_dir, '%s_%s_MEC_2018CV_vs_Valencia' % ( tag_neutrino_mec_hists, var ) ) )
429 
430  caption = open( os.path.join( plots_dir, '%s_%s_MEC_2018CV_vs_Valencia.txt' % ( tag_neutrino_mec_hists, var ) ), 'w' )
431  caption.write( text_caption + 'Shown are the 2018 CV MEC and Valencia MEC.' )
432  caption.close()
433 
434  legend3 = root_style.GetLegend()
435  legend3.SetTextSize( text_size )
436  legend3.SetTextFont( text_font )
437  #legend3.SetHeader( '%s CC MEC' % text_mec_neutrino )
438  legend3.AddEntry( hist_1D_tuned , text_nova_2p2h , 'L' )
439  legend3.AddEntry( hist_1D_minerva, text_minerva_2p2h, 'L' )
440  legend3.SetX1( xpos_legend_mec )
441  legend3.SetY2( ypos_legend_mec )
442  root_style.SetLegendPosY1( legend3 )
443 
444  canvas.Clear()
445 
446  hist_1D_tuned.SetMaximum( 1.8 * max( max_minerva, max_tuned ) )
447  hist_1D_tuned.Draw( 'hist' )
448  hist_1D_minerva.Draw( 'hist same' )
449  hist_1D_tuned.Draw( 'hist same' )
450  legend3.Draw()
451  if opts.sim : label_simulation.Draw()
452  label_beam.Draw()
453  label_numu_cc.Draw()
454  root_style.SavePlots( canvas, os.path.join( plots_dir, '%s_%s_MEC_2018CV_vs_MINERvA' % ( tag_neutrino_mec_hists, var ) ) )
455 
456  caption = open( os.path.join( plots_dir, '%s_%s_MEC_2018CV_vs_MINERvA.txt' % ( tag_neutrino_mec_hists, var ) ), 'w' )
457  caption.write( text_caption + 'Shown are the 2018 CV MEC and the MINERvA MEC tune.' )
458  caption.close()
459 
460  legend4 = root_style.GetLegend()
461  legend4.SetTextSize( text_size )
462  legend4.SetTextFont( text_font )
463  #legend4.SetHeader( '%s CC MEC' % text_mec_neutrino )
464  legend4.AddEntry( hist_1D_tuned , text_nova_2p2h , 'L' )
465  legend4.AddEntry( hist_1D_minerva , text_minerva_2p2h, 'L' )
466  legend4.AddEntry( hist_1D_valencia, text_valencia_mec, 'L' )
467  legend4.SetX1( xpos_legend_mec )
468  legend4.SetY2( ypos_legend_mec )
469  root_style.SetLegendPosY1( legend4 )
470 
471  canvas.Clear()
472 
473  hist_1D_tuned.SetMaximum( 2.0 * max( max_valencia, max_minerva, max_tuned ) )
474  hist_1D_tuned.Draw( 'hist' )
475  hist_1D_minerva.Draw( 'hist same' )
476  hist_1D_valencia.Draw( 'hist same' )
477  hist_1D_tuned.Draw( 'hist same' )
478  legend4.Draw()
479  if opts.sim : label_simulation.Draw()
480  label_beam.Draw()
481  label_numu_cc.Draw()
482  root_style.SavePlots( canvas, os.path.join( plots_dir, '%s_%s_MEC_2018CV_vs_All' % ( tag_neutrino_mec_hists, var ) ) )
483 
484  caption = open( os.path.join( plots_dir, '%s_%s_MEC_2018CV_vs_All.txt' % ( tag_neutrino_mec_hists, var ) ), 'w' )
485  caption.write( text_caption + 'Shown are the 2018 CV MEC, Valencia MEC, and the MINERvA MEC tune.' )
486  caption.close()
487 
488  legend5 = root_style.GetLegend()
489  legend5.SetTextSize( text_size )
490  legend5.SetTextFont( text_font )
491  #legend5.SetHeader( '%s CC MEC' % text_mec_neutrino )
492  legend5.AddEntry( hist_1D_tuned , text_nova_2p2h , 'L' )
493  legend5.AddEntry( hist_1D_shift_down, text_nova_2p2h + ' -1 #sigma', 'L' )
494  legend5.AddEntry( hist_1D_shift_up, text_nova_2p2h + ' +1 #sigma', 'L' )
495  #legend5.AddEntry( hist_1D_shift_down, '-1 #sigma', 'L' )
496  #legend5.AddEntry( hist_1D_shift_up, '+1 #sigma', 'L' )
497  legend5.SetX1( xpos_legend_mec )
498  legend5.SetY2( ypos_legend_mec )
499  root_style.SetLegendPosY1( legend5 )
500 
501  canvas.Clear()
502 
503  hist_1D_tuned.SetMaximum( 2.0 * max( max_shift_up, max_shift_down, max_tuned ) )
504  if var == 'TrueQ0' : hist_1D_tuned.GetXaxis().SetRangeUser( 0.0, 0.8 )
505  hist_1D_tuned.Draw( 'hist' )
506  hist_1D_shift_down.Draw( 'hist same' )
507  hist_1D_shift_up.Draw( 'hist same' )
508  hist_1D_tuned.Draw( 'hist same' )
509  legend5.Draw()
510  if opts.sim : label_simulation.Draw()
511  label_beam.Draw()
512  label_numu_cc.Draw()
513  root_style.SavePlots( canvas, os.path.join( plots_dir, '%s_%s_MEC_Shifts' % ( tag_neutrino_mec_hists, var ) ) )
514 
515  caption = open( os.path.join( plots_dir, '%s_%s_MEC_Shifts.txt' % ( tag_neutrino_mec_hists, var ) ), 'w' )
516  caption.write( text_caption + 'Shown are the 2018 CV MEC and 2018 MEC (q_0,|q|) +/-1 sigma shifts. The shifts affect both shape and normalization.' )
517  caption.close()
518 
519  legend6 = root_style.GetLegend()
520  legend6.SetTextSize( text_size )
521  legend6.SetTextFont( text_font )
522  legend6.AddEntry( hist_1D_non_mec , text_nova_tune, 'L' )
523  legend6.AddEntry( hist_1D_non_mec_shift_up , 'QE-Like Shifts', 'L' )
524  legend6.AddEntry( hist_1D_non_mec_shift_down, 'RES-Like Shifts', 'L' )
525  legend6.SetX1( xpos_legend_no_mec )
526  legend6.SetY2( ypos_legend_no_mec )
527  root_style.SetLegendPosY1( legend6 )
528 
529  canvas.Clear()
530 
531  ymax_non_mec = 2.0
532  if var == 'TrueQ0' : ymax_non_mec = 1.3
533 
534  hist_1D_non_mec.SetMaximum( ymax_non_mec * max( max_non_mec_shift_up, max_non_mec_shift_down, max_non_mec ) )
535  hist_1D_non_mec.GetXaxis().SetDecimals( kTRUE )
536  hist_1D_non_mec.Draw( 'hist' )
537  hist_1D_non_mec_shift_down.Draw( 'hist same' )
538  hist_1D_non_mec_shift_up.Draw( 'hist same' )
539  hist_1D_non_mec.Draw( 'hist same' )
540  legend6.Draw()
541  if opts.sim : label_simulation.Draw()
542  label_beam.Draw()
543  label_no_mec.Draw()
544  root_style.SavePlots( canvas, os.path.join( plots_dir, '%s_%s_NonMEC_Shifts' % ( opts.beam.upper(), var ) ) )
545 
546  caption = open( os.path.join( plots_dir, '%s_%s_NonMEC_Shifts.txt' % ( opts.beam.upper(), var ) ), 'w' )
547  caption.write( 'The %s near detector %s distribution for events that pass the 2018 numu selection and are not true numu/numubar CC MEC. ' % ( opts.beam.upper(), text_variable ) )
548  caption.write( 'Shown are the 2018 CV and the "QE-Like" and "RES-Like" shifts used to fit the 2018 MEC (q_0,|q|) +/-1 sigma shifts. ' )
549  caption.write( 'See the 2018 cross-section tuning tech note (DocDB 27755) for more details.' )
550  caption.close()
551 
552 
553  SetRangeQ0Q3( hist_nominal )
554  SetRangeQ0Q3( hist_tuned )
555 
556  event_scale_power_2d = 3
557  event_scale_2d = 1.0 / pow( 10.0, event_scale_power_2d )
558  text_events_2d = '10^{%d} Events' % event_scale_power_2d
559  scale_hist_2d = mc_pot_scale * event_scale_2d
560 
561  hist_nominal.Scale( scale_hist_2d )
562  hist_tuned .Scale( scale_hist_2d )
563 
564  hist_nominal.GetZaxis().SetTitle( text_events_2d )
565  hist_tuned .GetZaxis().SetTitle( text_events_2d )
566  hist_weights.GetZaxis().SetTitle( 'Weight' )
567 
568  canvas.SetLeftMargin(0.13)
569  canvas.SetRightMargin(0.17)
570  #gROOT.ForceStyle()
571 
572  text_weights_neutrino = '#nu'
573  if isRHC : text_weights_neutrino = '#bar{#nu}'
574 
575  text_weights_1 = "Weights Applied to"
576  text_weights_2 = "%s Empirical MEC" % text_weights_neutrino
577  text_empirical = "%s CC Empirical MEC" % text_mec_neutrino
578  text_nova_2p2h = "%s CC %s 2p2h" % ( text_mec_neutrino, text_nova )
579 
580  xpos_label_2d = canvas.GetLeftMargin() + 0.03 # text_size
581  ypos_label_weights_1 = ypos_label_top
582  ypos_label_weights_2 = ypos_label_top - text_size
583  ypos_label_nova_2p2h = ypos_label_top - text_size
584  ypos_label_empirical = ypos_label_top - text_size - 0.005
585  #ypos_label_mec_2d = ypos_label_top - text_size * 2 - 0.01
586  #if not isRHC : ypos_label_mec_2d -= 0.005
587 
588  label_weights_1 = root_style.GetLabel( text_weights_1, xpos_label_2d, ypos_label_weights_1, text_size, kWhite, text_font, 13, text_angle, setNDC )
589  label_weights_2 = root_style.GetLabel( text_weights_2, xpos_label_2d, ypos_label_weights_2, text_size, kWhite, text_font, 13, text_angle, setNDC )
590  label_empirical = root_style.GetLabel( text_empirical, xpos_label_2d, ypos_label_empirical, text_size, kBlack, text_font, 13, text_angle, setNDC )
591  label_nova_2p2h = root_style.GetLabel( text_nova_2p2h, xpos_label_2d, ypos_label_nova_2p2h, text_size, kBlack, text_font, 13, text_angle, setNDC )
592  label_beam.SetX( xpos_label_2d )
593  #label_mec .SetX( xpos_label_2d )
594  #label_mec .SetY( ypos_label_mec_2d )
595 
596  for i in range( 1, hist_weights.GetNbinsX() + 1 ) :
597  for j in range( 1, hist_weights.GetNbinsY() + 1 ) :
598  if hist_weights.GetBinContent( i, j ) < 0 :
599  hist_weights.SetBinContent( i, j, 0.0 )
600 
601  #hist_weights.GetZaxis().SetRangeUser( 0, 5.0 )
602  hist_weights.GetZaxis().SetRangeUser( -1.0e-6, 5.0 )
603 
605 
606  hist_weights.Draw( 'colz' )
607  label_weights_1.Draw()
608  label_weights_2.Draw()
609  if opts.sim : label_simulation.Draw()
610  root_style.SavePlots( canvas, os.path.join( plots_dir, '%s_TrueQ0_vs_TrueQ3_MEC_Weights' % tag_neutrino_mec_weights ) )
611 
612  gStyle.SetPalette( kBird )
613 
614  text_neutrino = 'neutrino'
615  if isRHC : text_neutrino = 'antineutrino'
616 
617  caption = open( os.path.join( plots_dir, '%s_TrueQ0_vs_TrueQ3_MEC_Weights.txt' % tag_neutrino_mec_weights ), 'w' )
618  caption.write( 'The 2018 %s MEC CV weights in true energy transfer, q_0, vs. true three-momentum transfer magnitude, |q|, bins. ' % text_neutrino )
619  caption.write( 'Non-zero weights in the kinematically forbidden region, q_0 > |q|, are an artifact of smoothing the weights.' )
620  caption.close()
621 
622  zaxis_tag = "_FixedZ"
623  zmax = 10
624  if isRHC : zmax = 3
625 
626  hist_nominal.Draw( 'colz' )
627  label_beam.Draw()
628  label_empirical.Draw()
629  if opts.sim : label_simulation.Draw()
630  root_style.SavePlots( canvas, os.path.join( plots_dir, tag_neutrino_mec_hists + "_TrueQ0_vs_TrueQ3_MEC_Empirical" ) )
631 
632  hist_nominal.GetZaxis().SetRangeUser( 0, zmax )
633  hist_nominal.Draw( 'colz' )
634  label_beam.Draw()
635  label_empirical.Draw()
636  if opts.sim : label_simulation.Draw()
637  root_style.SavePlots( canvas, os.path.join( plots_dir, tag_neutrino_mec_hists + "_TrueQ0_vs_TrueQ3_MEC_Empirical" + zaxis_tag ) )
638 
639  text_caption = 'The %s near detector true energy transfer, q_0, vs. true three-momentum transfer magnitude, |q|, ' % opts.beam.upper()
640  text_caption += 'distribution for true %s CC MEC events that pass the 2018 numu selection. ' % tag_neutrino_mec_hists.lower()
641 
642  caption = open( os.path.join( plots_dir, '%s_TrueQ0_vs_TrueQ3_MEC_Empirical.txt' % tag_neutrino_mec_hists ), 'w' )
643  caption.write( text_caption + 'Shown is the GENIE Empirical MEC.' )
644  caption.close()
645 
646  hist_tuned.Draw( 'colz' )
647  label_beam.Draw()
648  label_nova_2p2h.Draw()
649  if opts.sim : label_simulation.Draw()
650  root_style.SavePlots( canvas, os.path.join( plots_dir, tag_neutrino_mec_hists + "_TrueQ0_vs_TrueQ3_MEC_2018CV" ) )
651 
652  hist_tuned.GetZaxis().SetRangeUser( 0, zmax )
653  hist_tuned.Draw( 'colz' )
654  label_beam.Draw()
655  label_nova_2p2h.Draw()
656  if opts.sim : label_simulation.Draw()
657  root_style.SavePlots( canvas, os.path.join( plots_dir, tag_neutrino_mec_hists + "_TrueQ0_vs_TrueQ3_MEC_2018CV" + zaxis_tag ) )
658 
659  caption = open( os.path.join( plots_dir, '%s_TrueQ0_vs_TrueQ3_MEC_2018CV.txt' % tag_neutrino_mec_hists ), 'w' )
660  caption.write( text_caption + 'Shown is the 2018 CV MEC.' )
661  caption.close()
662 
void split(double tt, double *fr)
def SetHistLineAtt(hist, width, color)
Eigen::ArrayXd ProjectionX(const Eigen::MatrixXd &mat)
Helper for Unweighted.
correl_xv GetYaxis() -> SetDecimals()
constexpr T pow(T x)
Definition: pow.h:75
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
gargamelle SetTitle("Gargamelle #nu_{e} CC data")
correl_xv GetXaxis() -> SetDecimals()
procfile open("FD_BRL_v0.txt")
Eigen::ArrayXd ProjectionY(const Eigen::MatrixXd &mat)
Helper for WeightingVariable.
ratio_hxv GetZaxis() -> SetRangeUser(0, 2)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
simulatedPeak Scale(1/simulationNormalisationFactor)