diff --git a/limits/physicsModels/HHModel.py b/limits/physicsModels/HHModel.py index e81435d84..c2e35d603 100644 --- a/limits/physicsModels/HHModel.py +++ b/limits/physicsModels/HHModel.py @@ -8,13 +8,28 @@ # HHModel -> implements the interfaces to combine ################################### - from HiggsAnalysis.CombinedLimit.PhysicsModel import * +from HiggsAnalysis.CombinedLimit.PhysicsModel import SM_HIGG_DECAYS +from HiggsAnalysis.CombinedLimit.PhysicsModel import SM_HIGG_PROD +from HiggsAnalysis.CombinedLimit.SMHiggsBuilder import SMHiggsBuilder from sympy import * from numpy import matrix from numpy import linalg from sympy import Matrix +processSource_dict= { + "ttH":["tth"], + "ggH":["ggh","ggF"], + "qqH":["qqh","VBF"], + "VH":["vh"], + "ZH":["zh"], + "WH":["wh"], + "tHq":["thq","THQ"], + "tHW":["thw","THW"], + "ggHH":["gghh"], + "qqHH":["qqhh"], +} + class VBFHHSample: def __init__(self, val_CV, val_C2V, val_kl, val_xs, label): @@ -177,7 +192,10 @@ def calculatecoeffients(self): class HHModel(PhysicsModel): """ Models the HH production as linear sum of 6 components (VBF) and 3 components (GGF) """ - def __init__(self, ggf_sample_list, vbf_sample_list, name): + #scaleBR=True to scale for branching ratio dependence + #setBR1 to overwrite the single Higgs decay mode and one of the double Higgs decay modes + #setBR2 to overwrite the second decay mode of double Higgs + def __init__(self, ggf_sample_list, vbf_sample_list, name, setBR1="", setBR2=""): PhysicsModel.__init__(self) self.name = name @@ -188,7 +206,28 @@ def __init__(self, ggf_sample_list, vbf_sample_list, name): self.ggf_formula = GGFHHFormula(ggf_sample_list) self.vbf_formula = VBFHHFormula(vbf_sample_list) + self.setBR1 = setBR1 + self.setBR2 = setBR2 + self.scaleBR = True + + # include branching ratio uncertainties in the systematics + # if only one single HH channel is considered, it is suggested to + # include the specific BR uncertainty(s) directly in the datacard. + # This should be re-evaluated for a HH combination + self.doBRU = False + self.dump_inputs() + self.f_r_singleH_names=[] + + def setPhysicsOptions(self,physOptions): + for po in physOptions: + if po.startswith("BRU="): + self.doBRU = (po.replace("BRU=","") in [ "yes", "1", "Yes", "True", "true" ]) + if po.startswith("scaleBR="): + self.scaleBR = (po.replace("BRU=","") in [ "no", "0", "No", "False", "false" ]) + + print "BR uncertainties in partial widths: %s " % self.doBRU + def check_validity_ggf(self, ggf_sample_list): if len(ggf_sample_list) != 3: @@ -235,6 +274,7 @@ def doParametersOfInterest(self): self.modelBuilder.doVar("kt[1,-10,10]") self.modelBuilder.doSet("POI",POIs) + self.SMH = SMHiggsBuilder(self.modelBuilder) self.modelBuilder.out.var("r_gghh") .setConstant(True) self.modelBuilder.out.var("r_qqhh") .setConstant(True) @@ -243,11 +283,111 @@ def doParametersOfInterest(self): self.modelBuilder.out.var("kl") .setConstant(True) self.modelBuilder.out.var("kt") .setConstant(True) + #I need to build MH variables because the BR are tabulated as a function of MH + # the mass setting must be provided as input, i.e. '-m 125' + if self.modelBuilder.out.var("MH"): + self.modelBuilder.out.var("MH").setVal(self.options.mass) + self.modelBuilder.out.var("MH").setConstant(True) + else: + self.modelBuilder.doVar("MH[%g]" % self.options.mass) + self.create_scalings() def create_scalings(self): """ create the functions that scale the six components of vbf and the 3 components of ggf """ + ###################################################################### + #create Higgs BR scalings + for d in SM_HIGG_DECAYS + [ "hss" ]: + self.SMH.makeBR(d) + + # BR uncertainties + if self.doBRU: + self.SMH.makePartialWidthUncertainties() + else: + for d in SM_HIGG_DECAYS: + self.modelBuilder.factory_('HiggsDecayWidth_UncertaintyScaling_%s[1.0]' % d) + + # fix to have all BRs add up to unity + self.modelBuilder.factory_("sum::c7_SMBRs(%s)" % (",".join("SM_BR_"+X for X in "hzz hww htt hmm hcc hbb hss hgluglu hgg hzg".split()))) + self.modelBuilder.out.function("c7_SMBRs").Print("") + + # get VBF, tHq, tHW, ggZH cross section and resolved loops + self.SMH.makeScaling('qqH', CW='CV', CZ='CV') + self.SMH.makeScaling("tHq", CW='CV', Ctop="kt") + self.SMH.makeScaling("tHW", CW='CV', Ctop="kt") + self.SMH.makeScaling("ggZH",CZ='CV', Ctop="kt",Cb="1") + self.SMH.makeScaling('ggH', Cb='1', Ctop='kt', Cc="1") + self.SMH.makeScaling('hgluglu', Cb='1', Ctop='kt') + self.SMH.makeScaling('hgg', Cb='1', Ctop='kt', CW='CV', Ctau='1') + self.SMH.makeScaling('hzg', Cb='1', Ctop='kt', CW='CV', Ctau='1') + + + cGammap = {"hgg":0.49e-2,"hzz":0.83e-2,"hww":0.73e-2,"hgluglu":0.66e-2,"htt":0,"hbb":0,"hcc":0,"hmm":0} + + # First we need to create the terms that account for the self-coupling --> Just scale partial width first - https://arxiv.org/abs/1709.08649 Eq 22. + # probably a better way to code this since the partial width expressions are being repeated when we write the BR + for dec in cGammap.keys(): + valC1 = cGammap[dec] + self.modelBuilder.factory_('expr::kl_scalBR_%s("(@0-1)*%g",kl)' % (dec,valC1)) + + # next make the partial widths, also including the kappas -> we want to include the term from the normal kappas and the one from the self-coupling + self.modelBuilder.factory_('expr::kVktkl_Gscal_Z("(@0*@0+@3)*@1*@2", CV, SM_BR_hzz, HiggsDecayWidth_UncertaintyScaling_hzz, kl_scalBR_hzz)') + self.modelBuilder.factory_('expr::kVktkl_Gscal_W("(@0*@0+@3)*@1*@2", CV, SM_BR_hww, HiggsDecayWidth_UncertaintyScaling_hww, kl_scalBR_hww)') + self.modelBuilder.factory_('expr::kVktkl_Gscal_tau("(1+@4)*@0*@2 + (1+@5)*@1*@3", SM_BR_htt, SM_BR_hmm, HiggsDecayWidth_UncertaintyScaling_htt, HiggsDecayWidth_UncertaintyScaling_hmm,kl_scalBR_htt, kl_scalBR_hmm)') + self.modelBuilder.factory_('expr::kVktkl_Gscal_top("(1+@2)*@0*@1", SM_BR_hcc, HiggsDecayWidth_UncertaintyScaling_hcc, kl_scalBR_hcc)') + self.modelBuilder.factory_('expr::kVktkl_Gscal_bottom("(1+@3) * (@0*@2+@1)", SM_BR_hbb, SM_BR_hss, HiggsDecayWidth_UncertaintyScaling_hbb, kl_scalBR_hbb)') + self.modelBuilder.factory_('expr::kVktkl_Gscal_gluon(" (@0+@3) * @1 * @2", Scaling_hgluglu, SM_BR_hgluglu, HiggsDecayWidth_UncertaintyScaling_hgluglu, kl_scalBR_hgluglu)') + self.modelBuilder.factory_('expr::kVktkl_Gscal_gamma("(@0+@6)*@1*@4 + @2*@3*@5", Scaling_hgg, SM_BR_hgg, Scaling_hzg, SM_BR_hzg, HiggsDecayWidth_UncertaintyScaling_hgg, HiggsDecayWidth_UncertaintyScaling_hzg, kl_scalBR_hgg)') # no kl dependance on H->zg known yet ? + # fix to have all BRs add up to unity + self.modelBuilder.factory_("sum::kVktkl_SMBRs(%s)" % (",".join("SM_BR_"+X for X in "hzz hww htt hmm hcc hbb hss hgluglu hgg hzg".split()))) + self.modelBuilder.out.function("kVktkl_SMBRs").Print("") + + ## total witdh, normalized to the SM one (just the sum over the partial widths/SM total BR) + self.modelBuilder.factory_('expr::kVktkl_Gscal_tot("(@0+@1+@2+@3+@4+@5+@6)/@7", kVktkl_Gscal_Z, kVktkl_Gscal_W, kVktkl_Gscal_tau, kVktkl_Gscal_top, kVktkl_Gscal_bottom, kVktkl_Gscal_gluon, kVktkl_Gscal_gamma, kVktkl_SMBRs)') + + ## BRs, normalized to the SM ones: they scale as (partial/partial_SM) / (total/total_SM) + self.modelBuilder.factory_('expr::kVktkl_BRscal_hww("(@0*@0+@3)*@2/@1", CV, kVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_hww, kl_scalBR_hww)') + self.modelBuilder.factory_('expr::kVktkl_BRscal_hzz("(@0*@0+@3)*@2/@1", CV, kVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_hzz, kl_scalBR_hzz)') + self.modelBuilder.factory_('expr::kVktkl_BRscal_htt("(1+@2)*@1/@0", kVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_htt, kl_scalBR_htt)') + self.modelBuilder.factory_('expr::kVktkl_BRscal_hmm("(1+@2)*@1/@0", kVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_hmm, kl_scalBR_hmm)') + self.modelBuilder.factory_('expr::kVktkl_BRscal_hbb("(1+@2)*@1/@0", kVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_hbb, kl_scalBR_hbb)') + self.modelBuilder.factory_('expr::kVktkl_BRscal_hcc("(1+@2)*@1/@0", kVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_hcc, kl_scalBR_hcc)') + self.modelBuilder.factory_('expr::kVktkl_BRscal_hgg("(@0+@3)*@2/@1", Scaling_hgg, kVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_hgg,kl_scalBR_hgg)') + self.modelBuilder.factory_('expr::kVktkl_BRscal_hzg("@0*@2/@1", Scaling_hzg, kVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_hzg)') + self.modelBuilder.factory_('expr::kVktkl_BRscal_hgluglu("(@0+@3)*@2/@1", Scaling_hgluglu, kVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_hgluglu, kl_scalBR_hgluglu)') + + ###################################################################### + #create single Higgs production scalings + energy="13TeV" + cXSmap_13 = {"ggH":0.66e-2,"qqH":0.64e-2,"WH":1.03e-2,"ZH":1.19e-2,"ttH":3.51e-2,"VH":(0.5*(1.03e-2+1.19e-2))} + EWKmap_13 = {"ggH":1.049,"qqH":0.932,"WH":0.93,"ZH":0.947,"ttH":1.014,"VH":(0.5*(0.93+0.947))} + dZH = -1.536e-3 + + for production in SM_HIGG_PROD: + + if production in [ "ggZH", "tHq", "tHW"]: + self.f_r_singleH_names.append("Scaling_%s_%s"%(production,energy)) + + elif production in [ "ggH", "qqH" ]: + EWK = EWKmap_13[production] + self.modelBuilder.factory_("expr::kVktkl_XSscal_%s_%s(\"(@1+(@0-1)*%g/%g)/((1-(@0*@0-1)*%g))\",kl,Scaling_%s_%s)"\ + %(production,energy,cXSmap_13[production],EWK,dZH,production,energy)) + self.f_r_singleH_names.append("kVktkl_XSscal_%s_%s" %(production,energy) ) + + elif production in [ "ZH", "WH", "VH"]: + EWK = EWKmap_13[production] + self.modelBuilder.factory_("expr::kVktkl_XSscal_%s_%s(\"(@1*@1+(@0-1)*%g/%g)/((1-(@0*@0-1)*%g))\",kl,CV)"\ + %(production,energy,cXSmap_13[production],EWK,dZH)) + self.f_r_singleH_names.append("kVktkl_XSscal_%s_%s" %(production,energy) ) + + elif production == "ttH": + EWK = EWKmap_13[production] + self.modelBuilder.factory_("expr::kVktkl_XSscal_%s_%s(\"(@1*@1+(@0-1)*%g/%g)/((1-(@0*@0-1)*%g))\",kl,kt)"\ + %(production,energy,cXSmap_13[production],EWK,dZH)) + self.f_r_singleH_names.append("kVktkl_XSscal_%s_%s" %(production,energy) ) + + self.f_r_vbf_names = [] # the RooFormulae that scale the components (VBF) self.f_r_ggf_names = [] # the RooFormulae that scale the components (GGF) @@ -342,14 +482,56 @@ def pow_to_mul_string(expr): self.f_r_vbf_names.append(f_prod_name) #bookkeep the scaling that has been created def getYieldScale(self,bin,process): + + if not self.DC.isSignal[process]: return 1 + processSource = process.split("_")[0] + # It might happen that I use a different process naming than the default one + # In that case I will convert it to the default naming + if not processSource in processSource_dict: + for defaultname, equivalentnames in processSource_dict.items(): + for equivalentname in equivalentnames: + if processSource==equivalentname: + print "[WARNING]: process name \"%s\" is not the default name but I found it equivalent to \"%s\""%(processSource,defaultname) + print " --> if this is not the case fix the naming!" + processSource=defaultname + + + decaySource = [] + if self.setBR1!="": + decaySource.append(self.setBR1) + else: + decaySource.append(process.split("_")[-1]) + + #in case of HH I have to consider an additional Higgs decay + if processSource == "ggHH" or processSource == "qqHH": + if self.setBR2!="": + decaySource.append(self.setBR2) + else: + decaySource.append(process.split("_")[-2]) + + if not (processSource == "ggHH" or processSource == "qqHH"): + if(processSource == "bbH"): return 1 #neglect bbH dependence on kl + for XS_scaling_name in self.f_r_singleH_names: + if processSource in XS_scaling_name: + if not self.modelBuilder.out.function("kVktkl_BRscal_"+decaySource[0]): + raise RuntimeError, "Decay mode %s not supported" % decaySource[0] + if not self.scaleBR: + return XS_scaling_name + else: + BR_scaling_name = "kVktkl_BRscal_"+decaySource[0] + XSBR_scaling_name = "kVktkl_XSBRscal_%s_%s" % (processSource,decaySource[0]) + if not self.modelBuilder.out.function(XSBR_scaling_name): + self.modelBuilder.factory_("expr::%s(\"0.+@0*(@0>0.)*@1\",%s,%s)"%(XSBR_scaling_name,XS_scaling_name,BR_scaling_name))#I used a trick to avoid negative XS reweights + return XSBR_scaling_name + + #if I am here I have a double H process or an unsopported process ## my control to verify for a unique association between process <-> scaling function try: self.scalingMap except AttributeError: self.scalingMap = {} - if not self.DC.isSignal[process]: return 1 # match the process name in the datacard to the input sample of the calculation # this is the only point where the two things must be matched @@ -377,15 +559,42 @@ def getYieldScale(self,bin,process): if len(imatched_ggf) == 1: isample = imatched_ggf[0] self.scalingMap[process].append((isample, 'GGF')) - return self.f_r_ggf_names[isample] + XS_scaling_name = self.f_r_ggf_names[isample] + if not self.scaleBR: + return XS_scaling_name + else: + BR1_scaling_name = "kVktkl_BRscal_"+decaySource[0] + BR2_scaling_name = "kVktkl_BRscal_"+decaySource[1] + if not self.modelBuilder.out.function(BR1_scaling_name): + raise RuntimeError, "Decay mode %s not supported" %BR1_scaling_name + if not self.modelBuilder.out.function(BR2_scaling_name): + raise RuntimeError, "Decay mode %s not supported" %BR2_scaling_name + XSBR_scaling_name = "%s_%s_%s"%(XS_scaling_name,decaySource[0],decaySource[1]) + if not self.modelBuilder.out.function(XSBR_scaling_name): + self.modelBuilder.factory_('expr::%s("@0*@1*@2",%s,%s,%s)'%(XSBR_scaling_name,XS_scaling_name,BR1_scaling_name,BR2_scaling_name)) + return XSBR_scaling_name if len(imatched_vbf) == 1: isample = imatched_vbf[0] self.scalingMap[process].append((isample, 'VBF')) - return self.f_r_vbf_names[isample] + XS_scaling_name = self.f_r_vbf_names[isample] + if not self.scaleBR: + return XS_scaling_name + else: + BR1_scaling_name = "kVktkl_BRscal_"+decaySource[0] + BR2_scaling_name = "kVktkl_BRscal_"+decaySource[1] + if not self.modelBuilder.out.function(BR1_scaling_name): + raise RuntimeError, "Decay mode %s not supported" %BR1_scaling_name + if not self.modelBuilder.out.function(BR2_scaling_name): + raise RuntimeError, "Decay mode %s not supported" %BR2_scaling_name + XSBR_scaling_name = "%s_%s_%s"%(XS_scaling_name,decaySource[0],decaySource[1]) + if not self.modelBuilder.out.function(XSBR_scaling_name): + self.modelBuilder.factory_('expr::%s("@0*@1*@2",%s,%s,%s)'%(XSBR_scaling_name,XS_scaling_name,BR1_scaling_name,BR2_scaling_name)) + return XSBR_scaling_name raise RuntimeError('HHModel : fatal error in getYieldScale - this should never happen') + def done(self): ## this checks that a scaling has been attached to a unique process scalings = {} @@ -398,8 +607,8 @@ def done(self): for key, val in scalings.items(): if len(val) > 1: - print "[ERROR] : in HH model named", self.name, "there is a double assignment of a scaling : ", key, " ==> ", val - raise RuntimeError('HHModel : coudl not uniquely match the scaling to the process') + print "[WARNING] : in HH model named", self.name, "there is a double assignment of a scaling : ", key, " ==> ", val + #raise RuntimeError('HHModel : coudl not uniquely match the scaling to the process') ## now check that, if a VBF/GGF scaling exists, there are actually 6/3 samples in the card n_VBF = 0 @@ -413,10 +622,10 @@ def done(self): else: raise RuntimeError("HHModel : unrecognised type %s - should never happen" % i[0][1]) - if n_GGF > 0 and n_GGF != 3: + if n_GGF > 0 and n_GGF < 3: raise RuntimeError("HHModel : you did not pass all the 3 samples needed to build the GGF HH model") - if n_VBF > 0 and n_VBF != 6: + if n_VBF > 0 and n_VBF < 6: raise RuntimeError("HHModel : you did not pass all the 6 samples needed to build the VBF HH model") ######################################################## @@ -428,27 +637,30 @@ def done(self): ## NOTE 2 : the xs *must* correspond to the generator one # VBF : val_CV, val_C2V, val_kl + VBF_sample_list = [ VBFHHSample(1,1,1, val_xs = 0.00054/(0.3364), label = 'qqHH_CV_1_C2V_1_kl_1' ), VBFHHSample(1,2,1, val_xs = 0.00472/(0.3364), label = 'qqHH_CV_1_C2V_2_kl_1' ), VBFHHSample(1,1,2, val_xs = 0.00044/(0.3364), label = 'qqHH_CV_1_C2V_1_kl_2' ), VBFHHSample(1,1,0, val_xs = 0.00145/(0.3364), label = 'qqHH_CV_1_C2V_1_kl_0' ), VBFHHSample(0.5,1,1, val_xs = 0.00353/(0.3364), label = 'qqHH_CV_0p5_C2V_1_kl_1'), - VBFHHSample(1.5,1,1, val_xs = 0.02149/(0.3364), label = 'qqHH_CV_1p5_C2V_1_kl_1'), + VBFHHSample(1.5,1,1, val_xs = 0.02149/(0.3364), label = 'qqHH_CV_1p5_C2V_1_kl_1') ] # VBF : val_kl, val_kt GGF_sample_list = [ GGFHHSample(1,1, val_xs = 0.02675, label = 'ggHH_kl_1_kt_1' ), GGFHHSample(0,1, val_xs = 0.06007, label = 'ggHH_kl_0_kt_1' ), - GGFHHSample(5,1, val_xs = 0.07903, label = 'ggHH_kl_5_kt_1' ), - # GGFHHSample(2.45,1, val_xs = 0.01133, label = 'ggHH_kl_2p45_kt_1' ), + GGFHHSample(5,1, val_xs = 0.07903, label = 'ggHH_kl_5_kt_1' ) + #GGFHHSample(2.45,1, val_xs = 0.01133, label = 'ggHH_kl_2p45_kt_1' ) ] HHdefault = HHModel( ggf_sample_list = GGF_sample_list, vbf_sample_list = VBF_sample_list, - name = 'HHdefault' + name = 'HHdefault', + setBR1 = "",#setBR1 = "hgg", + setBR2 = "" #setBR2 = "hbb" ) # g = GGFHHFormula(GGF_sample_list) @@ -459,4 +671,4 @@ def done(self): # s1 : s.sample_list[0].val_xs, # s2 : s.sample_list[1].val_xs, # s3 : s.sample_list[2].val_xs, -# }) \ No newline at end of file +# })