make_mec_shifts_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 line_width = 4
14 top_margin = 0.08
15 bottom_margin = 0.12
16 left_margin = 0.15
17 right_margin = 0.05
18 
19 
20 def GetMaximum( hist ) :
21 
22  return hist.GetBinContent( hist.GetMaximumBin() )
23 
24 def PrepHist( hist, scale, ytitle, line_style, color ) :
25 
26  root_style.ApplyAxisStyleAll( hist )
27  hist.Scale( scale )
28  hist.GetYaxis().SetTitle( ytitle )
29  hist.SetLineWidth( line_width )
30  hist.SetLineStyle( line_style )
31  hist.SetLineColor( color )
32  hist.SetMarkerColor( color )
33 
34 if __name__ == "__main__":
35 
36  ###############################
37  # Parse Options
38  ###############################
39 
40  #default_outdir = "/nova/app/users/mislivec/scripts_xs_tuning_blessed_plots/plots_q0q3_shifts_rpa_fix/"
41  default_outdir = "./plots/mec_shifts/"
42 
43  parser = OptionParser()
44  parser.add_option( "--beam" , dest = "beam" , default = None , help = "Beam (FHC or RHC)" )
45  parser.add_option( "--outdir" , dest = "outdir" , default = default_outdir, help = "Output plots directory:\n%s" % default_outdir )
46  parser.add_option( "--minerva", dest = "minerva", action = "store_true" , default = False, help = "Plot MINERvA 2p2h" )
47  parser.add_option( "--prelim" , dest = "prelim" , action = "store_true" , default = False, help = "Add \"NOvA Preliminary\"" )
48 
49  if len(sys.argv) < 2 : parser.parse_args( "--help".split() )
50 
51  ( opts, args ) = parser.parse_args()
52 
53  if not opts.beam : sys.exit( "Beam (FHC or RHC) requred" )
54 
55  #indir = "/pnfs/nova/persistent/users/mislivec/xs_tuning_hists/rpa_fix/"
56  #indir = "."
57  indir = "/pnfs/nova/persistent/users/mislivec/xsec_tuning_paper/"
58 
59  hists_file_name_data = os.path.join( indir, "hists_2018_%s_data_full.root" % opts.beam )
60  hists_file_name_mc = os.path.join( indir, "hists_2018_%s_mc_full.root" % opts.beam )
61  if not os.path.exists( hists_file_name_data ) : sys.exit( "File does not exist:\n%s" % hists_file_name_data )
62  if not os.path.exists( hists_file_name_mc ) : sys.exit( "File does not exist:\n%s" % hists_file_name_mc )
63  print "Input Data histogram file:\n", hists_file_name_data
64  print "Input MC histogram file:\n", hists_file_name_mc
65 
66  hists_file_data = TFile( hists_file_name_data )
67  hists_file_mc = TFile( hists_file_name_mc )
68 
69  tree_pot_data = hists_file_data.Get( "pot" )
70  pot_data = 0.0
71  for i in range( tree_pot_data.GetEntries() ) :
72  tree_pot_data.GetEntry( i )
73  pot_data += tree_pot_data.pot_data
74 
75  tree_pot_mc = hists_file_mc.Get( "pot" )
76  pot_mc = 0.0
77  for i in range( tree_pot_mc.GetEntries() ) :
78  tree_pot_mc.GetEntry( i )
79  pot_mc += tree_pot_mc.pot_mc
80 
81  mc_pot_scale = pot_data / pot_mc
82 
83  print "pot_data = %.2e" % pot_data
84  print "pot_mc = %.2e" % pot_mc
85  print "mc_pot_scale = %.3f" % mc_pot_scale
86 
87  event_scale_power = 4
88  event_scale = 1.0 / pow( 10.0, event_scale_power )
89  text_events = "10^{%d} Events" % event_scale_power
90 
91  '''
92  for k in hists_file_mc.GetListOfKeys() :
93  name = k.GetName()
94  if "RecoEhadVis" in name : print name
95  #if "Shifts" in name : print name
96  '''
97 
98  #sys.exit()
99 
100  isRHC = opts.beam.upper() == "RHC"
101 
102  plots_dir = opts.outdir
103  if not os.path.isdir( plots_dir ) : os.makedirs( plots_dir )
104 
105  label_text_size = 0.04
106  text_font = 42
107  text_angle = 0.0
108  setNDC = True
109 
110  ypos_top = 1.0 - top_margin - 0.02
111 
112  text_beam = "Neutrino Beam"
113  if isRHC : text_beam = "Antineutrino Beam"
114  label_beam = root_style.GetLabel( text_beam, left_margin + label_text_size, ypos_top, label_text_size, kBlack, text_font, 13, text_angle, setNDC )
115 
116  text_nova = "NOvA" #"NO#nuA"
117  text_minerva = "MINERvA" #"MINER#nuA"
118 
119  label_prelim = root_style.GetLabel( text_nova + " Preliminary", 1.0 - right_margin, 1.0 - top_margin + 0.02, 0.06, kBlue, text_font, 31, text_angle, setNDC )
120 
121  for var in ( "RecoEhadVis", "RecoQ3" ) :
122 
123  hist_data = hists_file_data.Get( "%s_Data" % var )
124 
125  hist_cv_mec = hists_file_mc.Get( "%s_MEC" % var )
126  hist_cv_tot_mc = hists_file_mc.Get( "%s_TotMC" % var )
127  hist_cv_non_mec = hist_cv_tot_mc.Clone( "%s_NonMEC" % var )
128  hist_cv_non_mec.Add( hist_cv_mec, -1.0 )
129 
130  hist_neg_shift_tot_mc = hists_file_mc.Get( "%s_MECShapeDown_TotMC" % var )
131  hist_pos_shift_tot_mc = hists_file_mc.Get( "%s_MECShapeUp_TotMC" % var )
132  hist_minerva_tot_mc = hists_file_mc.Get( "%s_MEC_MINERvA" % var )
133  hist_minerva_tot_mc.Add( hist_cv_non_mec )
134 
135  PrepHist( hist_data, event_scale, text_events, 1, kBlack )
136 
137  mc_scale = event_scale * mc_pot_scale
138  PrepHist( hist_cv_tot_mc , mc_scale, text_events, 1, kBlack )
139  PrepHist( hist_cv_non_mec , mc_scale, text_events, 1, kGray )
140  PrepHist( hist_neg_shift_tot_mc, mc_scale, text_events, 2, kMagenta+1 )
141  PrepHist( hist_pos_shift_tot_mc, mc_scale, text_events, 2, kRed )
142  PrepHist( hist_minerva_tot_mc , mc_scale, text_events, 1, kAzure+1 )
143 
144  hist_cv_non_mec.SetFillStyle( 1001 )
145  hist_cv_non_mec.SetFillColor( kGray )
146  hist_cv_non_mec.SetLineColor( kGray+1 )
147 
148  legend = root_style.GetLegend()
149  legend.SetTextSize( 0.035 )
150  legend.SetTextFont( text_font )
151  legend.AddEntry( hist_data , text_nova + " ND Data" , "LPE" )
152  legend.AddEntry( hist_cv_tot_mc , text_nova + " 2018 #nu + #bar{#nu} Tune", "L" )
153  if opts.minerva :
154  legend.AddEntry( hist_minerva_tot_mc, text_minerva + " 2p2h" , "L" )
155  legend.AddEntry( hist_neg_shift_tot_mc, text_nova + " 2p2h -1#sigma" , "L" )
156  legend.AddEntry( hist_pos_shift_tot_mc, text_nova + " 2p2h +1#sigma" , "L" )
157  legend.AddEntry( hist_cv_non_mec , "Non-2p2h" , "F" )
158 
159  legend.SetX1( 0.5 )
160  legend.SetY2( ypos_top )
161  legend.SetY1( legend.GetY2() - ( legend.GetTextSize() + 0.005 ) * legend.GetNRows() )
162 
163  xmax = hist_data.GetXaxis().GetXmax()
164  xmin = 0.0
165  xndiv = hist_data.GetXaxis().GetNdivisions()
166  ymax_scale = 1.5
167  if var == "RecoEhadVis" :
168  if isRHC :
169  xmax = 0.4
170  xndiv = 504
171  ymax_scale = 1.3
172  else :
173  xmax = 0.6
174  xndiv = 506
175  ymax_scale = 1.3
176  elif var == "RecoQ3" :
177  if isRHC :
178  xmax = 1.4
179  xndiv = 407
180  ymax_scale = 1.8
181  else :
182  xmax = 1.6
183  xndiv = 408
184  ymax_scale = 1.8
185 
186  ymax = GetMaximum( hist_data )
187  hist_cv_non_mec.SetMaximum( ymax_scale * ymax )
188  hist_cv_non_mec.SetMinimum( 0.0 + 1.0e-3 )
189  hist_cv_non_mec.GetXaxis().SetRangeUser( xmin, xmax )
190  hist_cv_non_mec.GetXaxis().SetNdivisions( xndiv )
191 
192  hist_data.SetMarkerSize( 0.6 )
193  print "hist_data.GetMarkerSize() =", hist_data.GetMarkerSize()
194 
195  canvas_name = "Canvas_" + var
196  canvas = TCanvas( canvas_name, canvas_name, 300, 300 )
197  canvas.Divide( 1, 2, 0.0, 0.0 )
198 
199  canvas.cd(1)
200 
201  hist_cv_non_mec.Draw( "hist ][" )
202  if opts.minerva : hist_minerva_tot_mc.Draw( "hist ][ same" )
203  hist_neg_shift_tot_mc.Draw( "hist ][ same" )
204  hist_pos_shift_tot_mc.Draw( "hist ][ same" )
205  hist_cv_tot_mc.Draw( "hist ][ same" )
206  hist_data.Draw( "E same" )
207 
208  hist_data_ratio = hist_data.Clone( "%s_Data_Ratio" % var )
209  hist_pos_shift_tot_mc_ratio = hist_pos_shift_tot_mc.Clone( "%s_ShiftsQELike_TotMC_Ratio" % var )
210  hist_neg_shift_tot_mc_ratio = hist_neg_shift_tot_mc.Clone( "%s_ShiftsRESLike_TotMC_Ratio" % var )
211  hist_cv_tot_mc_ratio = hist_cv_tot_mc.Clone( "%s_TotMC_Ratio" % var )
212  hist_minerva_tot_mc_ratio = hist_minerva_tot_mc.Clone( "%s_TotMC_MINERvA_Ratio" % var )
213 
214  hist_data_ratio.Divide( hist_data )
215  hist_pos_shift_tot_mc_ratio.Divide( hist_data )
216  hist_neg_shift_tot_mc_ratio.Divide( hist_data )
217  hist_cv_tot_mc_ratio.Divide( hist_data )
218  hist_minerva_tot_mc_ratio.Divide( hist_data )
219 
220  hist_data_ratio.SetLineWidth( line_width )
221  hist_data_ratio.SetLineStyle( 7 )
222  hist_data_ratio.SetLineColor( kBlack )
223 
224  hist_data_ratio.GetXaxis().SetDecimals( kTRUE )
225  hist_data_ratio.GetXaxis().SetRangeUser( xmin, xmax )
226  if var == "RecoQ3" : hist_data_ratio.GetXaxis().SetTitle( "True #left|#vec{q}#right| (GeV/c)" )
227  hist_data_ratio.GetXaxis().SetNdivisions( xndiv )
228  hist_data_ratio.GetYaxis().SetRangeUser( 0.7 + 1.0e-3, 1.3 - 1.0e-3 )
229  hist_data_ratio.GetYaxis().SetTitle( "Sim. / Data" )
230 
231  canvas.cd( 2 )
232 
233  hist_data_ratio.Draw( "hist ][" )
234  if opts.minerva : hist_minerva_tot_mc_ratio.Draw( "hist C same" )
235  hist_neg_shift_tot_mc_ratio.Draw( "hist C same" )
236  hist_pos_shift_tot_mc_ratio.Draw( "hist C same" )
237  hist_cv_tot_mc_ratio.Draw( "hist C same" )
238 
239  canvas.SetTopMargin( 0.0 )
240  canvas.SetBottomMargin( 0.0 )
241  canvas.SetLeftMargin( 0.0 )
242  canvas.SetRightMargin( 0.0 )
243 
244  lower_fraction = 0.35
245  upper_fraction = 1.0 - lower_fraction
246 
247  title_size = 0.05
248  label_size = 0.035
249 
250  pad_upper = canvas.GetPad( 1 )
251  pad_upper.SetPad( 0, lower_fraction, 1, 1 - canvas.GetTopMargin() )
252  pad_upper.SetTopMargin( top_margin / upper_fraction )
253  pad_upper.SetBottomMargin( 0.0 )
254  pad_upper.SetLeftMargin( left_margin )
255  pad_upper.SetRightMargin( right_margin )
256  pad_upper.SetTicks( 1, 1 )
257 
258  for obj in pad_upper.GetListOfPrimitives() :
259  if obj.InheritsFrom( "TH1" ) :
260  obj.GetYaxis().SetTitleSize( title_size / upper_fraction )
261  obj.GetYaxis().SetTitleOffset( obj.GetYaxis().GetTitleOffset() * upper_fraction )
262  obj.GetYaxis().SetLabelSize( label_size / upper_fraction )
263 
264  pad_lower = canvas.GetPad( 2 )
265  pad_lower.SetPad( 0, 0.0, 1, lower_fraction )
266  pad_lower.SetTopMargin( 0.0 )
267  pad_lower.SetBottomMargin( bottom_margin / lower_fraction )
268  pad_lower.SetLeftMargin( left_margin )
269  pad_lower.SetRightMargin( right_margin )
270  pad_lower.SetTicks( 1, 1 )
271 
272  for obj in pad_lower.GetListOfPrimitives() :
273  if obj.InheritsFrom( "TH1" ) :
274  obj.GetYaxis().SetTitleSize( title_size / lower_fraction )
275  obj.GetYaxis().SetTitleOffset( obj.GetYaxis().GetTitleOffset() * lower_fraction )
276 
277  obj.GetYaxis().SetLabelSize( label_size / lower_fraction )
278  obj.GetXaxis().SetLabelSize( label_size / lower_fraction )
279  obj.GetXaxis().SetLabelOffset( 0.005 / lower_fraction )
280 
281  obj.GetXaxis().SetTitleSize( title_size / lower_fraction )
282  obj.GetXaxis().SetTitleOffset( 0.0 )
283 
284  pad_lower.Update()
285 
286  canvas.cd( 0 )
287 
288  if opts.prelim : label_prelim.Draw()
289  label_beam.Draw()
290  legend.Draw()
291 
292  plot_name = "%s_%s_MEC_Systs" % ( opts.beam.upper(), var )
293  if opts.minerva : plot_name += "_MINERvA"
294 
295  canvas.SaveAs( os.path.join( plots_dir, plot_name + ".pdf" ) )
296 
void split(double tt, double *fr)
def PrepHist(hist, scale, ytitle, line_style, color)
constexpr T pow(T x)
Definition: pow.h:75
gargamelle SetTitle("Gargamelle #nu_{e} CC data")