make_syst_table_plots.py
Go to the documentation of this file.
1 import os, sys, shlex, subprocess
2 from optparse import OptionParser
3 from ROOT import *
4 
5 gROOT.SetBatch( kTRUE )
6 
7 import root_style
8 
9 gROOT.SetBatch( kTRUE )
10 root_style.SetRootEnv()
11 
12 
13 systs_xsec = (
14 "CCQE z-exp EV shift #1",
15 "CCQE z-exp EV shift #2",
16 "CCQE z-exp EV shift #3",
17 "CCQE z-exp EV shift #4",
18 "MaCCRES",
19 "MvCCRES",
20 "MaNCRES",
21 "MvNCRES",
22 "ZNormCCQE",
23 "RPA shape: higher-Q^{2} enhancement (2020)",
24 "RPA shape: low-Q^{2} suppression (2020)",
25 "RES low-Q^2 suppression",
26 "DIS vnCC1pi",
27 "hN FSI mean free path",
28 "hN FSI fate fraction eigenvector #1",
29 "MEC E_{#nu} shape, neutrinos",
30 "MEC E_{#nu} shape, antineutrinos",
31 "MEC 2020 (q_{0}, |#vec{q}|) response, neutrinos",
32 "MEC 2020 (q_{0}, |#vec{q}|) response, antineutrinos",
33 "MEC initial state np fraction, neutrinos",
34 "MEC initial state np fraction, antineutrinos",
35 "Radiative corrections for #nu_{e}",
36 "Radiative corrections for #bar{#nu}_{e}",
37 "Second class currents",
38 "Genie PC 0",
39 "Genie PC 1",
40 "Genie PC 2",
41 "Genie PC 3",
42 "Genie PC 4",
43 "#nu_{#tau} Scale",
44 )
45 
46 systs_other = (
47 "Flux Component 00",
48 "Flux Component 01",
49 "Flux Component 02",
50 "Flux Component 03",
51 "Flux Component 04",
52 "AbsCalib",
53 "RelCalib",
54 "CalibShape",
55 "CalibDrift",
56 "Lightlevel",
57 "Cherenkov",
58 "Uncorr ND Mu Energy Scale 2020",
59 "Uncorr MuCat Mu Energy 2020",
60 "Neutron Pile-up 2020",
61 "Corr Mu Energy Scale 2020",
62 "Uncorr FD Mu Energy Scale 2020",
63 "Lepton Angle Syst ND XZ",
64 "Lepton Angle Syst ND YZ",
65 "Lepton Angle Syst FD XZ",
66 "Lepton Angle Syst FD YZ",
67 "Neutron visible energy systematic 2018",
68 "Acceptance ND to FD Kinematics Signal FHC 2020",
69 "Acceptance ND to FD Kinematics Signal RHC 2020",
70 "Michel Electrons Tagging Uncertainty",
71 )
72 
73 
74 def ReplaceSystName( syst_name ) :
75 
76  syst_name = syst_name.replace( "antineutrinos", "#bar{#nu}" )
77  syst_name = syst_name.replace( "neutrinos", "#nu" )
78  syst_name = syst_name.replace( "CCQE z-exp EV shift", "CCQE Z-Exp EV" )
79  syst_name = syst_name.replace( "ZNormCCQE", "CCQE Z-Exp Norm." )
80  syst_name = syst_name.replace( "RPA shape: higher-Q^{2} enhancement (2020)", "CCQE RPA High-Q^{2} Enh." )
81  syst_name = syst_name.replace( "RPA shape: low-Q^{2} suppression (2020)", "CCQE RPA Low-Q^{2} Supp." )
82  syst_name = syst_name.replace( "RES low-Q^2 suppression", "RES Low-Q^{2} Supp." )
83  syst_name = syst_name.replace( "hN FSI mean free path", "hN FSI MFP" )
84  syst_name = syst_name.replace( "hN FSI fate fraction eigenvector #1", "hN FSI Fate Frac. EV #1" )
85  syst_name = syst_name.replace( "MEC E_{#nu} shape", "MEC E_{#nu} Shape" )
86  syst_name = syst_name.replace( "MEC 2020 (q_{0}, |#vec{q}|) response", "MEC (q_{0},#left|#vec{q}#right|)" )
87  syst_name = syst_name.replace( "Radiative corrections for", "Rad. Corr." )
88  syst_name = syst_name.replace( "MEC initial state np fraction", "MEC np Frac." )
89  syst_name = syst_name.replace( "AbsCalib", "Abs. Calib." )
90  syst_name = syst_name.replace( "RelCalib", "Rel. Calib." )
91  syst_name = syst_name.replace( "CalibShape", "Calib. Shape" )
92  syst_name = syst_name.replace( "CalibDrift", "Calib. Drift" )
93  syst_name = syst_name.replace( "Lightlevel ", "Light Level" )
94  syst_name = syst_name.replace( "Neutron visible energy systematic 2018", "Neutron Vis. E" )
95  syst_name = syst_name.replace( "Uncorr ND Mu Energy Scale 2020", "E_{#mu} Scale ND" )
96  syst_name = syst_name.replace( "Uncorr FD Mu Energy Scale 2020", "E_{#mu} Scale FD" )
97  syst_name = syst_name.replace( "Corr Mu Energy Scale 2020", "E_{#mu} Scale Corr" )
98  syst_name = syst_name.replace( "Uncorr MuCat Mu Energy 2020", "E_{#mu} #mu Catcher" )
99  syst_name = syst_name.replace( "Neutron Pile-up 2020", "E_{#mu} Neutron Pile-up ND" )
100  syst_name = syst_name.replace( "Acceptance ND to FD Kinematics Signal FHC 2020", "Nue Signal Accept. FHC" )
101  syst_name = syst_name.replace( "Acceptance ND to FD Kinematics Signal RHC 2020", "Nue Signal Accept. RHC" )
102  syst_name = syst_name.replace( "Michel Electrons Tagging Uncertainty", "Michel Tagging" )
103  syst_name = syst_name.replace( "Flux Component 0", "Flux PC " )
104 
105  return syst_name
106 
107 def GetSystTable( fname, systs, addLepAng = False ) :
108 
109  table = []
110  for line in open( fname ) :
111 
112  if "Source of Uncertainty" in line : continue
113  if "hline" in line : continue
114 
115  line = line.strip()
116  if not "Statistical" in line : line = line[ 0 : line.find( "\\" ) ]
117  line = line.split( " & " )
118  syst_name = line[0]
119 
120  if systs == "xsec" and syst_name not in systs_xsec : continue
121  elif systs == "other" and syst_name not in systs_other : continue
122 
123  table.append( line )
124 
125  # Ugly hack to add zero-filled lepton angle uncertainty for no-pT extrap
126  if addLepAng and "Detector Response" in syst_name :
127  table.append( ( "Lepton Angle", "+0 / -0", "+0 / -0", "+0 / -0" ) )
128  if addLepAng and "Uncorr FD Mu Energy Scale 2020" in syst_name :
129  table.append( ( "Lepton Angle Syst ND XZ", "+0 / -0", "+0 / -0", "+0 / -0" ) )
130  table.append( ( "Lepton Angle Syst ND YZ", "+0 / -0", "+0 / -0", "+0 / -0" ) )
131  table.append( ( "Lepton Angle Syst FD XZ", "+0 / -0", "+0 / -0", "+0 / -0" ) )
132  table.append( ( "Lepton Angle Syst FD YZ", "+0 / -0", "+0 / -0", "+0 / -0" ) )
133 
134  return table
135 
136 
137 if __name__ == "__main__" :
138 
139  ###############################
140  # Parse Options
141  ###############################
142 
143  systs_options = ( "summary", "xsec", "other" )
144 
145  parser = OptionParser()
146  parser.add_option( "--systs" , dest = "systs" , type = "str" , default = None , help = "Systs to plot: " + ", ".join( systs_options ) )
147  parser.add_option( "--compare", dest = "compare", action = "store_true", default = False, help = "Plot tables for pT vs. no pT extrap comparison" )
148  parser.add_option( "--pt" , dest = "pt" , action = "store_true", default = False, help = "Plot tables for pT extrap" )
149 
150  if len(sys.argv) < 2 : parser.parse_args( "--help".split() )
151 
152  ( opts, args ) = parser.parse_args()
153 
154  if not opts.systs : sys.exit( "Specify systs" )
155  if opts.systs not in systs_options : sys.exit( "Wrong systs" )
156 
157  systs = opts.systs
158 
159  isSummary = systs == "summary"
160  isXSec = systs == "xsec"
161  isOther = systs == "other"
162 
163  table_type = "full"
164  if isSummary : table_type = "short"
165 
166  box_colors = ( kRed+1, kOrange+1, kGreen+1, kBlue )
167 
168  systs_tables = []
169  out_tag = ""
170  if opts.pt :
171  out_tag = "pt"
172  systs_tables.append( GetSystTable( "tables/syst_table_ptExtrap_%s.txt" % table_type, systs ) )
173  elif opts.compare :
174  out_tag = "pt_vs_no_pt"
175  systs_tables.append( GetSystTable( "tables/syst_table_%s.txt" % table_type, systs, True ) )
176  systs_tables.append( GetSystTable( "tables/syst_table_ptExtrap_%s.txt" % table_type, systs ) )
177  else :
178  sys.exit( "Please specify \"pt\" or \"compare\"" )
179 
180  plots_dir = "./syst_plots"
181  if not os.path.isdir( plots_dir ) : os.makedirs( plots_dir )
182 
183  ntables = len( systs_tables )
184  nsysts = len( systs_tables[0] )
185 
186  if ntables > len( box_colors ) : sys.exit( "Fewer colors than tables" )
187 
188  c = TCanvas()
189  canvas_width = 1000
190  canvas_height = 800
191  if isXSec :
192  canvas_height = 1200
193  c.SetCanvasSize( canvas_width, canvas_height )
194  c.SetLeftMargin( 0.3 )
195  c.SetRightMargin( 0.1 )
196  c.SetTopMargin( 0.1 )
197  c.SetGridy( 1 )
198 
199  osc_vars = ( "s23", "dCP", "m32" )
200  for var in osc_vars :
201 
202  var_index = osc_vars.index( var ) + 1
203 
204  xtitle = ""
205  xmax = 1.0
206  ndiv = 502
207  if var == "s23" :
208  xtitle = "sin^{2}#theta_{23}"
209  if isXSec :
210  xmax = 0.004
211  ndiv = 402
212  elif isOther :
213  xmax = 0.02
214  ndiv = 402
215  else :
216  xmax = 0.04
217  ndiv = 402
218  elif var == "dCP" :
219  xtitle = "#delta_{CP} / #pi"
220  if isXSec :
221  xmax = 0.05
222  ndiv = 502
223  elif isOther :
224  xmax = 0.04
225  ndiv = 402
226  else :
227  xmax = 0.3
228  ndiv = 302
229  elif var == "m32" :
230  xtitle = "#Delta m^{2}_{32} (#times 10^{-3} eV^{2})"
231  if isXSec :
232  xmax = 0.01
233  elif isOther :
234  xmax = 0.02
235  ndiv = 502
236  else :
237  xmax = 0.05
238  ndiv = 502
239 
240  h = TH2F( var, "", 10, -xmax, +xmax, nsysts, 0, nsysts )
241 
242  boxes = []
243 
244  for syst_idx in range( nsysts ) :
245 
246  ibin = nsysts - syst_idx
247 
248  syst_name = systs_tables[0][syst_idx][0]
249  syst_name = ReplaceSystName( syst_name )
250 
251  h.GetYaxis().SetBinLabel( ibin, syst_name )
252 
253  for table_idx in range( len( systs_tables ) ) :
254 
255  errors = systs_tables[table_idx][syst_idx][var_index].split( "/" )
256 
257  xlo = float( errors[1] )
258  xhi = float( errors[0] )
259  space = 0.25
260  dy = ( 1.0 - space ) / ntables
261  yhi = ibin - ( space / 2 ) - table_idx * dy
262  ylo = yhi - dy
263  boxes.append( TBox( xlo, ylo, xhi, yhi ) )
264  boxes[-1].SetFillColor( box_colors[ table_idx ] )
265  boxes[-1].SetLineWidth( 0 )
266 
267 
268  line = TLine()
269  line.SetLineWidth( 1 )
270  line.SetLineStyle( 3 )
271  line.SetLineColor( kBlack )
272 
273  root_style.ApplyAxisStyleAll( h )
274  h.GetXaxis().SetTitle( "Uncertainty in " + xtitle )
275  if isXSec : h.GetYaxis().SetLabelSize( 0.035 )
276  h.GetXaxis().SetNdivisions( ndiv )
277  h.GetXaxis().SetNoExponent( kTRUE )
278  h.Draw()
279  for b in boxes : b.Draw()
280  line.DrawLine( 0, 0, 0, nsysts )
281  for outType in ( "pdf", "C" ) :
282  c.SaveAs( os.path.join( plots_dir, "syst_table_%s_%s_%s.%s" % ( var, out_tag, systs, outType ) ) )
283 
void split(double tt, double *fr)
bin1_2sigma SetFillColor(3)
def GetSystTable(fname, systs, addLepAng=False)
gargamelle SetTitle("Gargamelle #nu_{e} CC data")
hmean SetLineWidth(2)
procfile open("FD_BRL_v0.txt")