2 from array
import array
5 from optparse
import OptionParser
10 gROOT.SetBatch( kTRUE )
11 root_style.SetRootEnv()
20 hist.GetXaxis().SetRangeUser( 0.0, 1.0 )
21 hist.GetYaxis().SetRangeUser( 0.0, 0.8 )
27 return hist.GetBinContent( hist.GetMaximumBin() )
33 n_color_contours = 999
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 )
46 gStyle.SetNumberContours( n_color_contours )
47 gStyle.SetPalette( n_color_contours, colors )
51 hist.GetXaxis().
SetTitle(
"True #left|#vec{q}#right| (GeV/c)" )
55 hist.GetXaxis().SetDecimals( kTRUE )
56 hist.GetYaxis().SetDecimals( kTRUE )
57 hist.GetZaxis().SetDecimals( kTRUE )
61 hist.SetLineWidth( width )
62 hist.SetLineColor( color )
64 if __name__ ==
"__main__":
71 default_indir =
"/pnfs/nova/persistent/users/mislivec/xsec_tuning_paper/" 72 default_outdir =
'./plots/true_q0q3' 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" )
83 if len(sys.argv) < 2 : parser.parse_args(
'--help'.
split() )
85 ( opts, args ) = parser.parse_args()
87 if not opts.beam : sys.exit(
'Beam (FHC or RHC) requred' )
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 )
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
100 hists_file_data = TFile( hists_file_name_data )
101 hists_file_mc = TFile( hists_file_name_mc )
103 tree_pot_data = hists_file_data.Get(
'pot' )
105 for i
in range( tree_pot_data.GetEntries() ) :
106 tree_pot_data.GetEntry( i )
107 pot_data += tree_pot_data.pot_data
109 tree_pot_mc = hists_file_mc.Get(
'pot' )
111 for i
in range( tree_pot_mc.GetEntries() ) :
112 tree_pot_mc.GetEntry( i )
113 pot_mc += tree_pot_mc.pot_mc
115 mc_pot_scale = pot_data / pot_mc
117 print 'pot_data = %.2e' % pot_data
118 print 'pot_mc = %.2e' % pot_mc
119 print 'mc_pot_scale = %.3f' % mc_pot_scale
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
126 for k
in hists_file_mc.GetListOfKeys() :
133 isRHC = opts.beam.upper() ==
'RHC' 135 tag_neutrino_mec_weights =
'Nu' 136 if isRHC : tag_neutrino_mec_weights =
'Nubar' 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 )
141 print 'Input MEC weights file:' 142 print weights_file_name
146 weights_file = TFile( weights_file_name )
147 for k
in weights_file.GetListOfKeys() :
print k.GetName()
150 plots_dir = default_outdir
151 if not os.path.isdir( plots_dir ) : os.makedirs( plots_dir )
153 tag_neutrino_mec_hists =
'Numu' 154 if isRHC : tag_neutrino_mec_hists =
'Numubar' 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() )
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 ) 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 )
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 )
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 ) 236 margin_left = gStyle.GetPadLeftMargin()
237 margin_right = gStyle.GetPadRightMargin()
238 margin_top = gStyle.GetPadTopMargin()
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
252 ypos_label_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
257 text_mec_neutrino =
'#nu_{#mu}' 258 if isRHC : text_mec_neutrino =
'#bar{#nu}_{#mu}' 260 text_beam =
"Neutrino Beam" 261 if isRHC : text_beam =
"Antineutrino Beam" 264 text_minerva =
"MINERvA" 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
275 label_beam = root_style.GetLabel( text_beam , xpos_label_mec , ypos_label_top , 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 )
283 for var
in (
'TrueQ0',
'TrueQ3' ) :
287 hist_1D_nominal = hist_nominal .
ProjectionY(
'%s_Empirical' % var )
289 hist_1D_valencia = hist_valencia .
ProjectionY(
'%s_Valencia' % 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 )
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 )
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 )
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 )
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 )
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 )
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 )
351 hist_1D_non_mec_shift_down.GetYaxis().
SetTitle( text_events_1d )
361 SetHistLineAtt( hist_1D_non_mec_shift_up , hist_line_width, kBlue )
362 SetHistLineAtt( hist_1D_non_mec_shift_down, hist_line_width, kRed )
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 )
374 legend1 = root_style.GetLegend()
375 legend1.SetTextSize( text_size )
376 legend1.SetTextFont( text_font )
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 )
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' )
391 if opts.sim : label_simulation.Draw()
394 root_style.SavePlots( canvas, os.path.join( plots_dir,
'%s_%s_MEC_2018CV_vs_Empirical' % ( tag_neutrino_mec_hists, var ) ) )
396 text_variable =
'blah' 398 text_variable =
'true energy transfer, q_0,' 400 text_variable =
'true three-momentum transfer magnitude, |q|,' 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() )
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.' )
408 legend2 = root_style.GetLegend()
409 legend2.SetTextSize( text_size )
410 legend2.SetTextFont( text_font )
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 )
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' )
425 if opts.sim : label_simulation.Draw()
428 root_style.SavePlots( canvas, os.path.join( plots_dir,
'%s_%s_MEC_2018CV_vs_Valencia' % ( tag_neutrino_mec_hists, var ) ) )
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.' )
434 legend3 = root_style.GetLegend()
435 legend3.SetTextSize( text_size )
436 legend3.SetTextFont( text_font )
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 )
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' )
451 if opts.sim : label_simulation.Draw()
454 root_style.SavePlots( canvas, os.path.join( plots_dir,
'%s_%s_MEC_2018CV_vs_MINERvA' % ( tag_neutrino_mec_hists, var ) ) )
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.' )
460 legend4 = root_style.GetLegend()
461 legend4.SetTextSize( text_size )
462 legend4.SetTextFont( text_font )
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 )
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' )
479 if opts.sim : label_simulation.Draw()
482 root_style.SavePlots( canvas, os.path.join( plots_dir,
'%s_%s_MEC_2018CV_vs_All' % ( tag_neutrino_mec_hists, var ) ) )
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.' )
488 legend5 = root_style.GetLegend()
489 legend5.SetTextSize( text_size )
490 legend5.SetTextFont( text_font )
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' )
497 legend5.SetX1( xpos_legend_mec )
498 legend5.SetY2( ypos_legend_mec )
499 root_style.SetLegendPosY1( legend5 )
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' )
510 if opts.sim : label_simulation.Draw()
513 root_style.SavePlots( canvas, os.path.join( plots_dir,
'%s_%s_MEC_Shifts' % ( tag_neutrino_mec_hists, var ) ) )
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.' )
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 )
532 if var ==
'TrueQ0' : ymax_non_mec = 1.3
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' )
541 if opts.sim : label_simulation.Draw()
544 root_style.SavePlots( canvas, os.path.join( plots_dir,
'%s_%s_NonMEC_Shifts' % ( opts.beam.upper(), var ) ) )
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.' )
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
561 hist_nominal.Scale( scale_hist_2d )
562 hist_tuned .
Scale( scale_hist_2d )
564 hist_nominal.GetZaxis().
SetTitle( text_events_2d )
566 hist_weights.GetZaxis().
SetTitle(
'Weight' )
568 canvas.SetLeftMargin(0.13)
569 canvas.SetRightMargin(0.17)
572 text_weights_neutrino =
'#nu' 573 if isRHC : text_weights_neutrino =
'#bar{#nu}' 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 )
580 xpos_label_2d = canvas.GetLeftMargin() + 0.03
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
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 )
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 )
602 hist_weights.GetZaxis().SetRangeUser( -1.0e-6, 5.0 )
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 ) )
612 gStyle.SetPalette( kBird )
614 text_neutrino =
'neutrino' 615 if isRHC : text_neutrino =
'antineutrino' 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.' )
622 zaxis_tag =
"_FixedZ" 626 hist_nominal.Draw(
'colz' )
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" ) )
632 hist_nominal.GetZaxis().SetRangeUser( 0, zmax )
633 hist_nominal.Draw(
'colz' )
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 ) )
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()
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.' )
646 hist_tuned.Draw(
'colz' )
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" ) )
652 hist_tuned.GetZaxis().SetRangeUser( 0, zmax )
653 hist_tuned.Draw(
'colz' )
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 ) )
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.' )
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()
def SetDecimalsAllAxes(hist)
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
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)
def SetMECWeightsPalette()
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
simulatedPeak Scale(1/simulationNormalisationFactor)