CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Public Attributes
ModelTools.ModelBuilder Class Reference
Inheritance diagram for ModelTools.ModelBuilder:
ModelTools.ModelBuilderBase ModelTools.CountingModelBuilder

Public Member Functions

def __init__
 
def doCombination
 
def doExpectedEvents
 
def doIndividualModels
 
def doModel
 
def doModelConfigs
 
def doNuisances
 
def doObservables
 
def isShapeSystematic
 
def setPhysics
 
- Public Member Functions inherited from ModelTools.ModelBuilderBase
def __init__
 
def doComment
 
def doObj
 
def doSet
 
def doVar
 
def factory_
 

Public Attributes

 DC
 
 physics
 
- Public Attributes inherited from ModelTools.ModelBuilderBase
 options
 
 out
 

Detailed Description

This class defines the actual methods to build a model

Definition at line 57 of file ModelTools.py.

Constructor & Destructor Documentation

def ModelTools.ModelBuilder.__init__ (   self,
  datacard,
  options 
)

Definition at line 59 of file ModelTools.py.

59 
60  def __init__(self,datacard,options):
61  ModelBuilderBase.__init__(self,options)
self.DC = datacard

Member Function Documentation

def ModelTools.ModelBuilder.doCombination (   self) 
create model_s and model_b pdfs

Definition at line 265 of file ModelTools.py.

Referenced by ModelTools.ModelBuilder.doModel().

265 
266  def doCombination(self):
267  """create model_s and model_b pdfs"""
raise RuntimeError, "Not implemented in ModelBuilder"
def ModelTools.ModelBuilder.doExpectedEvents (   self)

Definition at line 199 of file ModelTools.py.

References abs, ModelTools.ModelBuilderBase.doComment(), ModelTools.ModelBuilderBase.doObj(), ModelTools.ModelBuilderBase.doVar(), ModelTools.ModelBuilder.isShapeSystematic(), relativeConstraints.keys, and max().

Referenced by ModelTools.ModelBuilder.doModel().

199 
200  def doExpectedEvents(self):
201  self.doComment(" --- Expected events in each bin, for each process ----")
202  for b in self.DC.bins:
203  for p in self.DC.exp[b].keys(): # so that we get only self.DC.processes contributing to this bin
204  # if it's a zero background, write a zero and move on
205  if self.DC.exp[b][p] == 0:
206  self.doVar("n_exp_bin%s_proc_%s[%g]" % (b, p, self.DC.exp[b][p]))
207  continue
208  # get model-dependent scale factor
209  scale = self.physics.getYieldScale(b,p)
210  if scale == 0:
211  self.doVar("n_exp_bin%s_proc_%s[%g]" % (b, p, 0))
212  continue
213  # collect multiplicative corrections
214  nominal = self.DC.exp[b][p]
215  gamma = None; # gamma normalization (if present, DC.exp[b][p] is ignored)
216  factors = [] # RooAbsReal multiplicative factors (including gmN)
217  logNorms = [] # (kappa, RooAbsReal) lnN (or lnN)
218  alogNorms = [] # (kappaLo, kappaHi, RooAbsReal) asymm lnN
219  if scale == 1:
220  pass
221  elif type(scale) == str:
222  factors.append(scale)
223  else:
224  raise RuntimeError, "Physics model returned something which is neither a name, nor 0, nor 1."
225  for (n,nofloat,pdf,args,errline) in self.DC.systs:
226  if pdf == "param" or pdf.startswith("shape"): continue
227  if not errline[b].has_key(p): continue
228  if errline[b][p] == 0.0: continue
229  if pdf.startswith("shape") and pdf.endswith("?"): # might be a lnN in disguise
230  if not self.isShapeSystematic(b,p,n): pdf = "lnN"
231  if pdf == "lnN" and errline[b][p] == 1.0: continue
232  if pdf == "lnN" or pdf == "lnU":
233  if type(errline[b][p]) == list:
234  elow, ehigh = errline[b][p];
235  alogNorms.append((elow, ehigh, n))
236  else:
237  logNorms.append((errline[b][p], n))
238  elif pdf == "gmM":
239  factors.append(n)
240  elif pdf == "trG" or pdf == "unif":
241  myname = "n_exp_shift_bin%s_proc_%s_%s" % (b,p,n)
242  self.doObj(myname, ROOFIT_EXPR, "'1+%f*@0', %s" % (errline[b][p], n));
243  factors.append(myname)
244  elif pdf == "gmN":
245  factors.append(n)
246  if abs(errline[b][p] * args[0] - self.DC.exp[b][p]) > max(0.05 * max(self.DC.exp[b][p],1), errline[b][p]):
247  raise RuntimeError, "Values of N = %d, alpha = %g don't match with expected rate %g for systematics %s " % (
248  args[0], errline[b][p], self.DC.exp[b][p], n)
249  if gamma != None:
250  raise RuntimeError, "More than one gmN uncertainty for the same bin and process (second one is %s)" % n
251  gamma = n; nominal = errline[b][p];
252  else: raise RuntimeError, "Unsupported pdf %s" % pdf
253  # optimize constants
254  if len(factors) + len(logNorms) + len(alogNorms) == 0:
255  self.doVar("n_exp_bin%s_proc_%s[%g]" % (b, p, self.DC.exp[b][p]))
256  else:
257  #print "Process %s of bin %s depends on:\n\tlog-normals: %s\n\tasymm log-normals: %s\n\tother factors: %s\n" % (p,b,logNorms, alogNorms, factors)
258  procNorm = ROOT.ProcessNormalization("n_exp_bin%s_proc_%s" % (b,p), "", nominal)
259  for kappa, thetaName in logNorms: procNorm.addLogNormal(kappa, self.out.function(thetaName))
260  for kappaLo, kappaHi, thetaName in alogNorms: procNorm.addAsymmLogNormal(kappaLo, kappaHi, self.out.function(thetaName))
261  for factorName in factors: procNorm.addOtherFactor(self.out.function(factorName))
self.out._import(procNorm)
abs
#define abs(x)
Definition: mlp_lapack.h:159
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
def ModelTools.ModelBuilder.doIndividualModels (   self) 
create pdf_bin<X> and pdf_bin<X>_bonly for each bin

Definition at line 262 of file ModelTools.py.

Referenced by ModelTools.ModelBuilder.doModel().

262 
263  def doIndividualModels(self):
264  """create pdf_bin<X> and pdf_bin<X>_bonly for each bin"""
raise RuntimeError, "Not implemented in ModelBuilder"
def ModelTools.ModelBuilder.doModel (   self)

Definition at line 65 of file ModelTools.py.

References ModelTools.ModelBuilder.doCombination(), ModelTools.ModelBuilder.doExpectedEvents(), ModelTools.ModelBuilder.doIndividualModels(), ModelTools.ModelBuilder.doModelConfigs(), ModelTools.ModelBuilder.doNuisances(), and ModelTools.ModelBuilder.doObservables().

65 
66  def doModel(self):
67  self.doObservables()
68  self.physics.doParametersOfInterest()
69  self.physics.preProcessNuisances(self.DC.systs)
70  self.doNuisances()
71  self.doExpectedEvents()
72  self.doIndividualModels()
73  self.doCombination()
74  if self.options.bin:
75  self.doModelConfigs()
76  if self.options.verbose > 1: self.out.Print("V")
77  if self.options.verbose > 2:
78  print "Wrote GraphVizTree of model_s to ",self.options.out+".dot"
self.out.pdf("model_s").graphVizTree(self.options.out+".dot", "\\n")
def ModelTools.ModelBuilder.doModelConfigs (   self)

Definition at line 268 of file ModelTools.py.

References CgiWriter.out, ModelTools.ModelBuilderBase.out, CgiReader.out, ThePEG::HepMCConverter< HepMCEventT, Traits >::Vertex.out, PhysicsTools::Calibration::ProcMultiply.out, and BeamHaloAnalyzer.out.

Referenced by ModelTools.ModelBuilder.doModel().

268 
269  def doModelConfigs(self):
270  if not self.options.bin: raise RuntimeException
271  if self.options.out == None: raise RuntimeException
272  for nuis,warn in self.DC.flatParamNuisances.iteritems():
273  if self.out.var(nuis): self.out.var(nuis).setAttribute("flatParam")
274  elif warn: stderr.write("Missing variable %s declared as flatParam\n" % nuis)
275  mc_s = ROOT.RooStats.ModelConfig("ModelConfig", self.out)
276  mc_b = ROOT.RooStats.ModelConfig("ModelConfig_bonly", self.out)
277  for (l,mc) in [ ('s',mc_s), ('b',mc_b) ]:
278  mc.SetPdf(self.out.pdf("model_"+l))
279  mc.SetParametersOfInterest(self.out.set("POI"))
280  mc.SetObservables(self.out.set("observables"))
281  if len(self.DC.systs): mc.SetNuisanceParameters(self.out.set("nuisances"))
282  if self.out.set("globalObservables"): mc.SetGlobalObservables(self.out.set("globalObservables"))
283  if self.options.verbose > 1: mc.Print("V")
284  self.out._import(mc, mc.GetName())
self.out.writeToFile(self.options.out)
def ModelTools.ModelBuilder.doNuisances (   self)

Definition at line 82 of file ModelTools.py.

References ModelTools.ModelBuilderBase.doComment(), ModelTools.ModelBuilderBase.doObj(), ModelTools.ModelBuilderBase.doSet(), ModelTools.ModelBuilderBase.doVar(), join(), max(), split, mathSSE.sqrt(), and makeHLTPrescaleTable.values.

Referenced by ModelTools.ModelBuilder.doModel().

82 
83  def doNuisances(self):
84  if len(self.DC.systs) == 0: return
85  self.doComment(" ----- nuisances -----")
86  globalobs = []
87  for (n,nofloat,pdf,args,errline) in self.DC.systs:
88  if pdf == "lnN" or pdf.startswith("shape"):
89  r = "-4,4" if pdf == "shape" else "-7,7"
90  self.doObj("%s_Pdf" % n, "Gaussian", "%s[%s], %s_In[0,%s], 1" % (n,r,n,r));
91  globalobs.append("%s_In" % n)
92  if self.options.bin:
93  self.out.var("%s_In" % n).setConstant(True)
94  elif pdf == "gmM":
95  val = 0;
96  for c in errline.values(): #list channels
97  for v in c.values(): # list effects in each channel
98  if v != 0:
99  if val != 0 and v != val:
100  raise RuntimeError, "Error: line %s contains two different uncertainties %g, %g, which is not supported for gmM" % (n,v,val)
101  val = v;
102  if val == 0: raise RuntimeError, "Error: line %s contains all zeroes"
103  theta = val*val; kappa = 1/theta
104  self.doObj("%s_Pdf" % n, "Gamma", "%s[1,%f,%f], %s_In[%g,%g,%g], %s_scaling[%g], 0" % (n, max(0.01,1-5*val), 1+5*val, n, kappa, 1, 2*kappa+4, n, theta))
105  globalobs.append("%s_In" % n)
106  if self.options.bin: self.out.var("%s_In" % n).setConstant(True)
107  elif pdf == "gmN":
108  if False:
109  # old version, that creates a poisson with a very large range
110  self.doObj("%s_Pdf" % n, "Poisson", "%s_In[%d,0,%d], %s[0,%d], 1" % (n,args[0],2*args[0]+5,n,2*args[0]+5))
111  else:
112  # new version, that creates a poisson with a narrower range (but still +/- 7 sigmas)
113  #print "Searching for bounds for",n,"poisson with obs",args[0]
114  minExp = args[0]+1 if args[0] > 0 else 0;
115  while (ROOT.TMath.Poisson(args[0], minExp) > 1e-12) and minExp > 0:
116  #print "Poisson(%d, minExp = %f) = %g > 1e-12" % (args[0], minExp, ROOT.TMath.Poisson(args[0], minExp))
117  minExp *= 0.8;
118  maxExp = args[0]+1;
119  while (ROOT.TMath.Poisson(args[0], maxExp) > 1e-12):
120  #print "Poisson(%d, maxExp = %f) = %g > 1e-12" % (args[0], maxExp, ROOT.TMath.Poisson(args[0], maxExp))
121  maxExp *= 1.2;
122  minObs = args[0];
123  while minObs > 0 and (ROOT.TMath.Poisson(minObs, args[0]+1) > 1e-12):
124  #print "Poisson(minObs = %d, %f) = %g > 1e-12" % (minObs, args[0]+1, ROOT.TMath.Poisson(minObs, args[0]+1))
125  minObs -= (sqrt(args[0]) if args[0] > 10 else 1);
126  maxObs = args[0]+2;
127  while (ROOT.TMath.Poisson(maxObs, args[0]+1) > 1e-12):
128  #print "Poisson(maxObs = %d, %f) = %g > 1e-12" % (maxObs, args[0]+1, ROOT.TMath.Poisson(maxObs, args[0]+1))
129  maxObs += (sqrt(args[0]) if args[0] > 10 else 2);
130  self.doObj("%s_Pdf" % n, "Poisson", "%s_In[%d,%f,%f], %s[%f,%f,%f], 1" % (n,args[0],minObs,maxObs,n,args[0]+1,minExp,maxExp))
131  globalobs.append("%s_In" % n)
132  if self.options.bin: self.out.var("%s_In" % n).setConstant(True)
133  elif pdf == "trG":
134  trG_min = -7; trG_max = +7;
135  for b in errline.keys():
136  for v in errline[b].values():
137  if v > 0 and 1.0 + trG_min * v < 0: trG_min = -1.0/v;
138  if v < 0 and 1.0 + trG_max * v < 0: trG_max = -1.0/v;
139  r = "%f,%f" % (trG_min,trG_max);
140  self.doObj("%s_Pdf" % n, "Gaussian", "%s[0,%s], %s_In[0,%s], 1" % (n,r,n,r));
141  globalobs.append("%s_In" % n)
142  if self.options.bin:
143  self.out.var("%s_In" % n).setConstant(True)
144  elif pdf == "lnU":
145  self.doObj("%s_Pdf" % n, "Uniform", "%s[-1,1]" % n);
146  elif pdf == "unif":
147  self.doObj("%s_Pdf" % n, "Uniform", "%s[%f,%f]" % (n,args[0],args[1]))
148  elif pdf == "param":
149  mean = float(args[0])
150  if "/" in args[1]:
151  sigmaL,sigmaR = args[1].split("/")
152  if sigmaL[0] != "-" or sigmaR[0] != "+": raise RuntimeError, "Asymmetric parameter uncertainties should be entered as -x/+y"
153  sigmaL = sigmaL[1:]; sigmaR = sigmaR[1:]
154  if len(args) == 3: # mean, sigma, range
155  self.doVar("%s%s" % (n,args[2]))
156  else:
157  sigma = float(args[1])
158  if self.out.var(n):
159  self.out.var(n).setConstant(False)
160  self.out.var(n).setRange(mean-4*float(sigmaL),mean+4*float(sigmaR))
161  else:
162  self.doVar("%s[%g,%g]" % (n, mean-4*float(sigmaL), mean+4*float(sigmaR)))
163  self.out.var(n).setVal(mean)
164  self.doObj("%s_Pdf" % n, "BifurGauss", "%s, %s_In[%s,%g,%g], %s, %s" % (n, n, args[0], self.out.var(n).getMin(), self.out.var(n).getMax(), sigmaL, sigmaR))
165  self.out.var("%s_In" % n).setConstant(True)
166  else:
167  if len(args) == 3: # mean, sigma, range
168  self.doVar("%s%s" % (n,args[2]))
169  else:
170  sigma = float(args[1])
171  if self.out.var(n):
172  self.out.var(n).setConstant(False)
173  self.out.var(n).setRange(mean-4*sigma, mean+4*sigma)
174  else:
175  self.doVar("%s[%g,%g]" % (n, mean-4*sigma, mean+4*sigma))
176  self.out.var(n).setVal(mean)
177  self.doObj("%s_Pdf" % n, "Gaussian", "%s, %s_In[%s,%g,%g], %s" % (n, n, args[0], self.out.var(n).getMin(), self.out.var(n).getMax(), args[1]))
178  self.out.var("%s_In" % n).setConstant(True)
179  globalobs.append("%s_In" % n)
180  else: raise RuntimeError, "Unsupported pdf %s" % pdf
181  if nofloat:
182  self.out.var("%s" % n).setAttribute("globalConstrained",True)
183  if self.options.bin:
184  nuisPdfs = ROOT.RooArgList()
185  nuisVars = ROOT.RooArgSet()
186  for (n,nf,p,a,e) in self.DC.systs:
187  nuisVars.add(self.out.var(n))
188  nuisPdfs.add(self.out.pdf(n+"_Pdf"))
189  self.out.defineSet("nuisances", nuisVars)
190  self.out.nuisPdf = ROOT.RooProdPdf("nuisancePdf", "nuisancePdf", nuisPdfs)
191  self.out._import(self.out.nuisPdf)
192  self.out.nuisPdfs = nuisPdfs
193  gobsVars = ROOT.RooArgSet()
194  for g in globalobs: gobsVars.add(self.out.var(g))
195  self.out.defineSet("globalObservables", gobsVars)
196  else: # doesn't work for too many nuisances :-(
197  self.doSet("nuisances", ",".join(["%s" % n for (n,nf,p,a,e) in self.DC.systs]))
198  self.doObj("nuisancePdf", "PROD", ",".join(["%s_Pdf" % n for (n,nf,p,a,e) in self.DC.systs]))
self.doSet("globalObservables", ",".join(globalobs))
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
join
static std::string join(char **cmd)
Definition: RemoteFile.cc:18
split
double split
Definition: MVATrainer.cc:139
def ModelTools.ModelBuilder.doObservables (   self) 
create pdf_bin<X> and pdf_bin<X>_bonly for each bin

Definition at line 79 of file ModelTools.py.

Referenced by ModelTools.ModelBuilder.doModel().

79 
80  def doObservables(self):
81  """create pdf_bin<X> and pdf_bin<X>_bonly for each bin"""
raise RuntimeError, "Not implemented in ModelBuilder"
def ModelTools.ModelBuilder.isShapeSystematic (   self,
  channel,
  process,
  syst 
)

Definition at line 285 of file ModelTools.py.

Referenced by ModelTools.ModelBuilder.doExpectedEvents(), ShapeTools.ShapeBuilder.getExtraNorm(), and ShapeTools.ShapeBuilder.getPdf().

285 
286  def isShapeSystematic(self,channel,process,syst):
287  return False
def ModelTools.ModelBuilder.setPhysics (   self,
  physicsModel 
)

Definition at line 62 of file ModelTools.py.

62 
63  def setPhysics(self,physicsModel):
64  self.physics = physicsModel
self.physics.setModelBuilder(self)

Member Data Documentation

ModelTools.ModelBuilder.DC

Definition at line 61 of file ModelTools.py.

ModelTools.ModelBuilder.physics

Definition at line 63 of file ModelTools.py.

Generated for CMSSW Reference Manual by doxygen 1.8.5