Logo ROOT   6.12/07
Reference Guide
makeQuickModel.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 
3 #
4 # A pyROOT script that allows one to
5 # make quick measuremenst.
6 #
7 # This is a command-line script that
8 # takes in signal and background values,
9 # as well as potentially uncertainties on those
10 # values, and returns a fitted signal value
11 # and errors
12 
13 
14 def main():
15  """ Create a simple model and run statistical tests
16 
17  This script can be used to make simple statistical using histfactory.
18  It takes values for signal, background, and data as input, and
19  can optionally take uncertainties on signal or background.
20  The model is created and saved to an output ROOT file, and
21  the model can be fit if requested.
22 
23  """
24 
25  # Let's parse the input
26  # Define the command line options of the script:
27  import optparse
28  desc = " ".join(main.__doc__.split())
29 
30  vers = "0.1"
31  parser = optparse.OptionParser( description = desc, version = vers, usage = "%prog [options]" )
32 
33  parser.add_option( "-s", "--signal", dest = "signal",
34  action = "store", type = "float", default=None,
35  help = "Expected Signal" )
36 
37  parser.add_option( "-b", "--background", dest = "background",
38  action = "store", type = "float", default=None,
39  help = "Expected Background" )
40 
41  parser.add_option( "-d", "--data", dest = "data",
42  action = "store", type = "float", default=None,
43  help = "Measured data" )
44 
45  parser.add_option( "--signal-uncertainty", dest = "signal_uncertainty",
46  action = "store", type = "float", default=None,
47  help = "Uncertainty on the signal rate, as a fraction. --signal-uncertainty=.05 means a 5% uncertainty." )
48 
49  parser.add_option( "--background-uncertainty", dest = "background_uncertainty",
50  action = "store", type = "float", default=None,
51  help = "Uncertainty on the background rate, as a fraction, not a percentage. --background-uncertainty=.05 means a 5% uncertainty." )
52 
53  parser.add_option( "--output-prefix", dest = "output_prefix",
54  action = "store", type = "string", default="Measurement",
55  help = "Prefix for output files when using export. Can include directories (ie 'MyDirectory/MyPrefix')" )
56 
57  parser.add_option( "-e", "--export", dest = "export",
58  action = "store_true", default=False,
59  help = "Make output plots, graphs, and save the workspace." )
60 
61  # Parse the command line options:
62  ( options, unknown ) = parser.parse_args()
63 
64  # Make a log
65  # Set the format of the log messages:
66  FORMAT = 'Py:%(name)-25s %(levelname)-8s %(message)s'
67  import logging
68  logging.basicConfig( format = FORMAT )
69  # Create the logger object:
70  logger = logging.getLogger( "makeQuickMeasurement" )
71  # Set the following to DEBUG when debugging the scripts:
72  logger.setLevel( logging.INFO )
73 
74  # So a small sanity check:
75  if len( unknown ):
76  logger.warning( "Options(s) not recognised: [" + ",".join( unknown ) + "]" )
77 
78  # Ensure that all necessary input has been supplied
79  if options.signal == None:
80  logger.error( "You have to define a value for expacted signal (use --signal)" )
81  return 255
82 
83  if options.background == None:
84  logger.error( "You have to define a value for expacted background (use --background)" )
85  return 255
86 
87  if options.data == None:
88  logger.error( "You have to define a value for measured data (use --data)" )
89  return 255
90 
91 
92  # Okay, if all input is acceptable, we simply pass
93  # it to the MakeSimpleMeasurement function, which
94  # does the real work.
95 
96  MakeSimpleMeasurement( signal_val=options.signal, background_val=options.background, data_val=options.data,
97  signal_uncertainty=options.signal_uncertainty, background_uncertainty=options.background_uncertainty,
98  Export=options.export, output_prefix=options.output_prefix)
99  return
100 
101 
102 def MakeSimpleMeasurement( signal_val, background_val, data_val, signal_uncertainty=None, background_uncertainty=None,
103  Export=False, output_prefix="Measurement"):
104  """ Make a simple measurement using HistFactory
105 
106  Take in simple values for signal, background data,
107  and potentially uncertainty on signal and background
108 
109  """
110 
111  try:
112  import ROOT
113  except:
114  print "It seems that pyROOT isn't properly configured"
115  return
116 
117  # Create and name a measurement
118  # Set the Parameter of interest, and set several
119  # other parameters to be constant
120  meas = ROOT.RooStats.HistFactory.Measurement("meas", "meas")
121  meas.SetOutputFilePrefix( output_prefix )
122  meas.SetPOI( "SigXsecOverSM" )
123 
124  # We don't include Lumi here,
125  # but since HistFactory gives it to
126  # us for free, we set it constant
127  # The values are just dummies
128 
129  meas.SetLumi( 1.0 )
130  meas.SetLumiRelErr( 0.10 )
131  meas.AddConstantParam("Lumi")
132 
133  # We set ExportOnly to false. This parameter
134  # defines what happens when we run MakeMeasurementAndModelFast
135  # If we DO run that function, we also want it to export.
136  meas.SetExportOnly( False )
137 
138  # Create a channel and set
139  # the measured value of data
140  # (no extenal hist necessar for cut-and-count)
141  chan = ROOT.RooStats.HistFactory.Channel( "channel" )
142  chan.SetData( data_val )
143 
144  # Create the signal sample and set it's value
145  signal = ROOT.RooStats.HistFactory.Sample( "signal" )
146  signal.SetNormalizeByTheory( False )
147  signal.SetValue( signal_val )
148  #signal.SetValue( 10 )
149 
150  # Add the parmaeter of interest and a systematic
151  # Try to make intelligent choice of upper bound
152  import math
153  upper_bound = 3*math.ceil( (data_val - background_val) / signal_val )
154  upper_bound = max(upper_bound, 3)
155  signal.AddNormFactor( "SigXsecOverSM", 1, 0, upper_bound )
156 
157  # If we have a signal uncertainty, add it too
158  if signal_uncertainty != None:
159  uncertainty_up = 1.0 + signal_uncertainty
160  uncertainty_down = 1.0 - signal_uncertainty
161  signal.AddOverallSys( "signal_uncertainty", uncertainty_down, uncertainty_up )
162 
163  # Finally, add this sample to the channel
164  chan.AddSample( signal )
165 
166  # Create a background sample
167  background = ROOT.RooStats.HistFactory.Sample( "background" )
168  background.SetNormalizeByTheory( False )
169  background.SetValue( background_val )
170 
171  # If we have a background uncertainty, add it too
172  if background_uncertainty != None:
173  uncertainty_up = 1.0 + background_uncertainty
174  uncertainty_down = 1.0 - background_uncertainty
175  background.AddOverallSys( "background_uncertainty", uncertainty_down, uncertainty_up )
176 
177  # Finally, add this sample to the channel
178  chan.AddSample( background )
179 
180  # Add this channel to the measurement
181  # There is only this one channel, after all
182  meas.AddChannel( chan )
183 
184  # Now, do the measurement
185  if Export:
186  workspace = ROOT.RooStats.HistFactory.MakeModelAndMeasurementFast( meas )
187  return workspace
188 
189  else:
190  factory = ROOT.RooStats.HistFactory.HistoToWorkspaceFactoryFast()
191  workspace = factory.MakeCombinedModel( meas )
192  #workspace = ROOT.RooStats.HistFactory.HistoToWorkspaceFactoryFast.MakeCombinedModel( meas )
193  ROOT.RooStats.HistFactory.FitModel( workspace )
194 
195  # At this point, we are done
196  return
197 
198 
199 if __name__ == "__main__":
200  main()
def MakeSimpleMeasurement(signal_val, background_val, data_val, signal_uncertainty=None, background_uncertainty=None, Export=False, output_prefix="Measurement")