|
|
|
Public Member Functions | |
| def | __init__ |
| def | doCombination |
| def | doExpectedEvents |
| def | doIndividualModels |
| def | doModel |
| def | doModelConfigs |
| def | doNuisances |
| def | doObservables |
| def | isShapeSystematic |
| def | setPhysics |
Public Attributes | |
| DC | |
| physics | |
This class defines the actual methods to build a model
Definition at line 57 of file ModelTools.py.
| def ModelTools::ModelBuilder::__init__ | ( | self, | |
| datacard, | |||
| options | |||
| ) |
Reimplemented in ModelTools::CountingModelBuilder.
Definition at line 59 of file ModelTools.py.
| def ModelTools::ModelBuilder::doCombination | ( | self | ) |
create model_s and model_b pdfs
Reimplemented in ModelTools::CountingModelBuilder.
Definition at line 265 of file ModelTools.py.
| def ModelTools::ModelBuilder::doExpectedEvents | ( | self | ) |
Definition at line 199 of file ModelTools.py.
00200 : 00201 self.doComment(" --- Expected events in each bin, for each process ----") 00202 for b in self.DC.bins: 00203 for p in self.DC.exp[b].keys(): # so that we get only self.DC.processes contributing to this bin 00204 # if it's a zero background, write a zero and move on 00205 if self.DC.exp[b][p] == 0: 00206 self.doVar("n_exp_bin%s_proc_%s[%g]" % (b, p, self.DC.exp[b][p])) 00207 continue 00208 # get model-dependent scale factor 00209 scale = self.physics.getYieldScale(b,p) 00210 if scale == 0: 00211 self.doVar("n_exp_bin%s_proc_%s[%g]" % (b, p, 0)) 00212 continue 00213 # collect multiplicative corrections 00214 nominal = self.DC.exp[b][p] 00215 gamma = None; # gamma normalization (if present, DC.exp[b][p] is ignored) 00216 factors = [] # RooAbsReal multiplicative factors (including gmN) 00217 logNorms = [] # (kappa, RooAbsReal) lnN (or lnN) 00218 alogNorms = [] # (kappaLo, kappaHi, RooAbsReal) asymm lnN 00219 if scale == 1: 00220 pass 00221 elif type(scale) == str: 00222 factors.append(scale) 00223 else: 00224 raise RuntimeError, "Physics model returned something which is neither a name, nor 0, nor 1." 00225 for (n,nofloat,pdf,args,errline) in self.DC.systs: 00226 if pdf == "param" or pdf.startswith("shape"): continue 00227 if not errline[b].has_key(p): continue 00228 if errline[b][p] == 0.0: continue 00229 if pdf.startswith("shape") and pdf.endswith("?"): # might be a lnN in disguise 00230 if not self.isShapeSystematic(b,p,n): pdf = "lnN" 00231 if pdf == "lnN" and errline[b][p] == 1.0: continue 00232 if pdf == "lnN" or pdf == "lnU": 00233 if type(errline[b][p]) == list: 00234 elow, ehigh = errline[b][p]; 00235 alogNorms.append((elow, ehigh, n)) 00236 else: 00237 logNorms.append((errline[b][p], n)) 00238 elif pdf == "gmM": 00239 factors.append(n) 00240 elif pdf == "trG" or pdf == "unif": 00241 myname = "n_exp_shift_bin%s_proc_%s_%s" % (b,p,n) 00242 self.doObj(myname, ROOFIT_EXPR, "'1+%f*@0', %s" % (errline[b][p], n)); 00243 factors.append(myname) 00244 elif pdf == "gmN": 00245 factors.append(n) 00246 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]): 00247 raise RuntimeError, "Values of N = %d, alpha = %g don't match with expected rate %g for systematics %s " % ( 00248 args[0], errline[b][p], self.DC.exp[b][p], n) 00249 if gamma != None: 00250 raise RuntimeError, "More than one gmN uncertainty for the same bin and process (second one is %s)" % n 00251 gamma = n; nominal = errline[b][p]; 00252 else: raise RuntimeError, "Unsupported pdf %s" % pdf 00253 # optimize constants 00254 if len(factors) + len(logNorms) + len(alogNorms) == 0: 00255 self.doVar("n_exp_bin%s_proc_%s[%g]" % (b, p, self.DC.exp[b][p])) 00256 else: 00257 #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) 00258 procNorm = ROOT.ProcessNormalization("n_exp_bin%s_proc_%s" % (b,p), "", nominal) 00259 for kappa, thetaName in logNorms: procNorm.addLogNormal(kappa, self.out.function(thetaName)) 00260 for kappaLo, kappaHi, thetaName in alogNorms: procNorm.addAsymmLogNormal(kappaLo, kappaHi, self.out.function(thetaName)) 00261 for factorName in factors: procNorm.addOtherFactor(self.out.function(factorName)) self.out._import(procNorm)
| def ModelTools::ModelBuilder::doIndividualModels | ( | self | ) |
create pdf_bin<X> and pdf_bin<X>_bonly for each bin
Reimplemented in ModelTools::CountingModelBuilder.
Definition at line 262 of file ModelTools.py.
| def ModelTools::ModelBuilder::doModel | ( | self | ) |
Definition at line 65 of file ModelTools.py.
00066 : 00067 self.doObservables() 00068 self.physics.doParametersOfInterest() 00069 self.physics.preProcessNuisances(self.DC.systs) 00070 self.doNuisances() 00071 self.doExpectedEvents() 00072 self.doIndividualModels() 00073 self.doCombination() 00074 if self.options.bin: 00075 self.doModelConfigs() 00076 if self.options.verbose > 1: self.out.Print("V") 00077 if self.options.verbose > 2: 00078 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.
00269 : 00270 if not self.options.bin: raise RuntimeException 00271 if self.options.out == None: raise RuntimeException 00272 for nuis,warn in self.DC.flatParamNuisances.iteritems(): 00273 if self.out.var(nuis): self.out.var(nuis).setAttribute("flatParam") 00274 elif warn: stderr.write("Missing variable %s declared as flatParam\n" % nuis) 00275 mc_s = ROOT.RooStats.ModelConfig("ModelConfig", self.out) 00276 mc_b = ROOT.RooStats.ModelConfig("ModelConfig_bonly", self.out) 00277 for (l,mc) in [ ('s',mc_s), ('b',mc_b) ]: 00278 mc.SetPdf(self.out.pdf("model_"+l)) 00279 mc.SetParametersOfInterest(self.out.set("POI")) 00280 mc.SetObservables(self.out.set("observables")) 00281 if len(self.DC.systs): mc.SetNuisanceParameters(self.out.set("nuisances")) 00282 if self.out.set("globalObservables"): mc.SetGlobalObservables(self.out.set("globalObservables")) 00283 if self.options.verbose > 1: mc.Print("V") 00284 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.
00083 : 00084 if len(self.DC.systs) == 0: return 00085 self.doComment(" ----- nuisances -----") 00086 globalobs = [] 00087 for (n,nofloat,pdf,args,errline) in self.DC.systs: 00088 if pdf == "lnN" or pdf.startswith("shape"): 00089 r = "-4,4" if pdf == "shape" else "-7,7" 00090 self.doObj("%s_Pdf" % n, "Gaussian", "%s[%s], %s_In[0,%s], 1" % (n,r,n,r)); 00091 globalobs.append("%s_In" % n) 00092 if self.options.bin: 00093 self.out.var("%s_In" % n).setConstant(True) 00094 elif pdf == "gmM": 00095 val = 0; 00096 for c in errline.values(): #list channels 00097 for v in c.values(): # list effects in each channel 00098 if v != 0: 00099 if val != 0 and v != val: 00100 raise RuntimeError, "Error: line %s contains two different uncertainties %g, %g, which is not supported for gmM" % (n,v,val) 00101 val = v; 00102 if val == 0: raise RuntimeError, "Error: line %s contains all zeroes" 00103 theta = val*val; kappa = 1/theta 00104 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)) 00105 globalobs.append("%s_In" % n) 00106 if self.options.bin: self.out.var("%s_In" % n).setConstant(True) 00107 elif pdf == "gmN": 00108 if False: 00109 # old version, that creates a poisson with a very large range 00110 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)) 00111 else: 00112 # new version, that creates a poisson with a narrower range (but still +/- 7 sigmas) 00113 #print "Searching for bounds for",n,"poisson with obs",args[0] 00114 minExp = args[0]+1 if args[0] > 0 else 0; 00115 while (ROOT.TMath.Poisson(args[0], minExp) > 1e-12) and minExp > 0: 00116 #print "Poisson(%d, minExp = %f) = %g > 1e-12" % (args[0], minExp, ROOT.TMath.Poisson(args[0], minExp)) 00117 minExp *= 0.8; 00118 maxExp = args[0]+1; 00119 while (ROOT.TMath.Poisson(args[0], maxExp) > 1e-12): 00120 #print "Poisson(%d, maxExp = %f) = %g > 1e-12" % (args[0], maxExp, ROOT.TMath.Poisson(args[0], maxExp)) 00121 maxExp *= 1.2; 00122 minObs = args[0]; 00123 while minObs > 0 and (ROOT.TMath.Poisson(minObs, args[0]+1) > 1e-12): 00124 #print "Poisson(minObs = %d, %f) = %g > 1e-12" % (minObs, args[0]+1, ROOT.TMath.Poisson(minObs, args[0]+1)) 00125 minObs -= (sqrt(args[0]) if args[0] > 10 else 1); 00126 maxObs = args[0]+2; 00127 while (ROOT.TMath.Poisson(maxObs, args[0]+1) > 1e-12): 00128 #print "Poisson(maxObs = %d, %f) = %g > 1e-12" % (maxObs, args[0]+1, ROOT.TMath.Poisson(maxObs, args[0]+1)) 00129 maxObs += (sqrt(args[0]) if args[0] > 10 else 2); 00130 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)) 00131 globalobs.append("%s_In" % n) 00132 if self.options.bin: self.out.var("%s_In" % n).setConstant(True) 00133 elif pdf == "trG": 00134 trG_min = -7; trG_max = +7; 00135 for b in errline.keys(): 00136 for v in errline[b].values(): 00137 if v > 0 and 1.0 + trG_min * v < 0: trG_min = -1.0/v; 00138 if v < 0 and 1.0 + trG_max * v < 0: trG_max = -1.0/v; 00139 r = "%f,%f" % (trG_min,trG_max); 00140 self.doObj("%s_Pdf" % n, "Gaussian", "%s[0,%s], %s_In[0,%s], 1" % (n,r,n,r)); 00141 globalobs.append("%s_In" % n) 00142 if self.options.bin: 00143 self.out.var("%s_In" % n).setConstant(True) 00144 elif pdf == "lnU": 00145 self.doObj("%s_Pdf" % n, "Uniform", "%s[-1,1]" % n); 00146 elif pdf == "unif": 00147 self.doObj("%s_Pdf" % n, "Uniform", "%s[%f,%f]" % (n,args[0],args[1])) 00148 elif pdf == "param": 00149 mean = float(args[0]) 00150 if "/" in args[1]: 00151 sigmaL,sigmaR = args[1].split("/") 00152 if sigmaL[0] != "-" or sigmaR[0] != "+": raise RuntimeError, "Asymmetric parameter uncertainties should be entered as -x/+y" 00153 sigmaL = sigmaL[1:]; sigmaR = sigmaR[1:] 00154 if len(args) == 3: # mean, sigma, range 00155 self.doVar("%s%s" % (n,args[2])) 00156 else: 00157 sigma = float(args[1]) 00158 if self.out.var(n): 00159 self.out.var(n).setConstant(False) 00160 self.out.var(n).setRange(mean-4*float(sigmaL),mean+4*float(sigmaR)) 00161 else: 00162 self.doVar("%s[%g,%g]" % (n, mean-4*float(sigmaL), mean+4*float(sigmaR))) 00163 self.out.var(n).setVal(mean) 00164 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)) 00165 self.out.var("%s_In" % n).setConstant(True) 00166 else: 00167 if len(args) == 3: # mean, sigma, range 00168 self.doVar("%s%s" % (n,args[2])) 00169 else: 00170 sigma = float(args[1]) 00171 if self.out.var(n): 00172 self.out.var(n).setConstant(False) 00173 self.out.var(n).setRange(mean-4*sigma, mean+4*sigma) 00174 else: 00175 self.doVar("%s[%g,%g]" % (n, mean-4*sigma, mean+4*sigma)) 00176 self.out.var(n).setVal(mean) 00177 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])) 00178 self.out.var("%s_In" % n).setConstant(True) 00179 globalobs.append("%s_In" % n) 00180 else: raise RuntimeError, "Unsupported pdf %s" % pdf 00181 if nofloat: 00182 self.out.var("%s" % n).setAttribute("globalConstrained",True) 00183 if self.options.bin: 00184 nuisPdfs = ROOT.RooArgList() 00185 nuisVars = ROOT.RooArgSet() 00186 for (n,nf,p,a,e) in self.DC.systs: 00187 nuisVars.add(self.out.var(n)) 00188 nuisPdfs.add(self.out.pdf(n+"_Pdf")) 00189 self.out.defineSet("nuisances", nuisVars) 00190 self.out.nuisPdf = ROOT.RooProdPdf("nuisancePdf", "nuisancePdf", nuisPdfs) 00191 self.out._import(self.out.nuisPdf) 00192 self.out.nuisPdfs = nuisPdfs 00193 gobsVars = ROOT.RooArgSet() 00194 for g in globalobs: gobsVars.add(self.out.var(g)) 00195 self.out.defineSet("globalObservables", gobsVars) 00196 else: # doesn't work for too many nuisances :-( 00197 self.doSet("nuisances", ",".join(["%s" % n for (n,nf,p,a,e) in self.DC.systs])) 00198 self.doObj("nuisancePdf", "PROD", ",".join(["%s_Pdf" % n for (n,nf,p,a,e) in self.DC.systs])) self.doSet("globalObservables", ",".join(globalobs))
| def ModelTools::ModelBuilder::doObservables | ( | self | ) |
create pdf_bin<X> and pdf_bin<X>_bonly for each bin
Reimplemented in ModelTools::CountingModelBuilder.
Definition at line 79 of file ModelTools.py.
| def ModelTools::ModelBuilder::isShapeSystematic | ( | self, | |
| channel, | |||
| process, | |||
| syst | |||
| ) |
Definition at line 285 of file ModelTools.py.
| def ModelTools::ModelBuilder::setPhysics | ( | self, | |
| physicsModel | |||
| ) |
Definition at line 62 of file ModelTools.py.
Definition at line 59 of file ModelTools.py.
Definition at line 62 of file ModelTools.py.
1.7.3