Single h br scaling

• scale single-Higgs production processes as a function of CV, kl, kt with the option --PO 'doHscaling=True'
• with 'doHscaling=True' option, also the Higgs processes of background (positive process number) are scaled
• scale Higgs branching ratios as a function of CV, kl, kt with the option --PO 'doBRscaling=True'

For a correct production and decay parsing, the process names must be provided with the following syntax:
<production>_<whatever>_<decay>
For intermediate decays, e.g. H->ZZ->llll, the syntax must be <hdecay1><sub-decays>
Some example:

qqHH_CV_1_C2V_1_kl_0_2016_hbbhgg      
ggHH_kl_1_kt_1_hbbhzz4l     
ttH_hzz4l      
ggH_hww2l2nu      

dhi/models/HBRscaler.py

+
+# reorganization of https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/blob/102x/python/TrilinearCouplingModels.py
+# for HH combination
+
+from HiggsAnalysis.CombinedLimit.PhysicsModel import SM_HIGG_DECAYS,SM_HIGG_PROD
+from HiggsAnalysis.CombinedLimit.SMHiggsBuilder import SMHiggsBuilder
+
+# single Higgs production scalings coefficients 
+# WH and ZH coeff are very similar --> build VH coeff as an average btw the two 
+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
+
+class HBRscaler():
+    """ Produce single H and BR scalings for anomalous couplings, and produce XS*BR scalings for both H and HH"""
+
+    def __init__(self, modelBuilder, doBRscaling, doHscaling):
+        self.modelBuilder = modelBuilder
+        #use SMHiggsBuilder to build single H XS and BR scalings
+        self.SMH = SMHiggsBuilder(self.modelBuilder)
+        self.doBRscaling = doBRscaling
+        self.doHscaling = doHscaling
+        self.f_BR_scalings=[]
+        self.f_H_scalings=[]
+        self.nowarning_processes=[] #processes inside this list won't print warnings in case the scaling is not found
+
+        if self.doBRscaling:
+            self.BuildBRscalings()
+        if self.doHscaling:
+            self.BuildHscalings()
+
+    def BuildBRscalings(self):
+        
+        for d in SM_HIGG_DECAYS + [ "hss" ]:
+            self.SMH.makeBR(d)
+
+        # FIXME: to check how to deal with BR uncertainties --> for now keep them frozen
+        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("")        
+
+        # define resolved loops
+        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')
+        
+        # BR scaling vs kl (https://arxiv.org/abs/1709.08649 Eq 22)
+        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}
+        for dec in cGammap.keys(): 
+            valC1 = cGammap[dec]
+            self.modelBuilder.factory_('expr::kl_scalBR_%s("(@0-1)*%g",kl)' % (dec,valC1))
+
+        #partial widths as a function of kl, kt, and CV
+        self.modelBuilder.factory_('expr::CVktkl_Gscal_Z("(@0*@0+@3)*@1*@2", CV, SM_BR_hzz, HiggsDecayWidth_UncertaintyScaling_hzz, kl_scalBR_hzz)')
+        self.modelBuilder.factory_('expr::CVktkl_Gscal_W("(@0*@0+@3)*@1*@2", CV, SM_BR_hww, HiggsDecayWidth_UncertaintyScaling_hww, kl_scalBR_hww)')
+        self.modelBuilder.factory_('expr::CVktkl_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::CVktkl_Gscal_top("(1+@2)*@0*@1", SM_BR_hcc, HiggsDecayWidth_UncertaintyScaling_hcc, kl_scalBR_hcc)')
+        self.modelBuilder.factory_('expr::CVktkl_Gscal_bottom("(1+@3) * (@0*@2+@1)", SM_BR_hbb, SM_BR_hss, HiggsDecayWidth_UncertaintyScaling_hbb, kl_scalBR_hbb)')
+        self.modelBuilder.factory_('expr::CVktkl_Gscal_gluon("  (@0+@3)  * @1 * @2", Scaling_hgluglu, SM_BR_hgluglu, HiggsDecayWidth_UncertaintyScaling_hgluglu, kl_scalBR_hgluglu)')
+        self.modelBuilder.factory_('expr::CVktkl_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::CVktkl_SMBRs(%s)" %  (",".join("SM_BR_"+X for X in "hzz hww htt hmm hcc hbb hss hgluglu hgg hzg".split())))
+        #self.modelBuilder.out.function("CVktkl_SMBRs").Print("")        
+
+        ## total witdh, normalized to the SM one (just the sum over the partial widths/SM total BR)
+        self.modelBuilder.factory_('expr::CVktkl_Gscal_tot("(@0+@1+@2+@3+@4+@5+@6)/@7", CVktkl_Gscal_Z, CVktkl_Gscal_W, CVktkl_Gscal_tau, CVktkl_Gscal_top, CVktkl_Gscal_bottom, CVktkl_Gscal_gluon, CVktkl_Gscal_gamma, CVktkl_SMBRs)') 
+
+        ## BRs, normalized to the SM ones: they scale as (partial/partial_SM) / (total/total_SM) 
+        self.modelBuilder.factory_('expr::CVktkl_BRscal_hww("(@0*@0+@3)*@2/@1", CV, CVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_hww, kl_scalBR_hww)')
+        self.modelBuilder.factory_('expr::CVktkl_BRscal_hzz("(@0*@0+@3)*@2/@1", CV, CVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_hzz, kl_scalBR_hzz)')
+        self.modelBuilder.factory_('expr::CVktkl_BRscal_htt("(1+@2)*@1/@0", CVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_htt, kl_scalBR_htt)')
+        self.modelBuilder.factory_('expr::CVktkl_BRscal_hmm("(1+@2)*@1/@0", CVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_hmm, kl_scalBR_hmm)')
+        self.modelBuilder.factory_('expr::CVktkl_BRscal_hbb("(1+@2)*@1/@0", CVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_hbb, kl_scalBR_hbb)')
+        self.modelBuilder.factory_('expr::CVktkl_BRscal_hcc("(1+@2)*@1/@0", CVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_hcc, kl_scalBR_hcc)')
+        self.modelBuilder.factory_('expr::CVktkl_BRscal_hgg("(@0+@3)*@2/@1", Scaling_hgg, CVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_hgg,kl_scalBR_hgg)')
+        self.modelBuilder.factory_('expr::CVktkl_BRscal_hzg("@0*@2/@1", Scaling_hzg, CVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_hzg)')
+        self.modelBuilder.factory_('expr::CVktkl_BRscal_hgluglu("(@0+@3)*@2/@1", Scaling_hgluglu, CVktkl_Gscal_tot, HiggsDecayWidth_UncertaintyScaling_hgluglu, kl_scalBR_hgluglu)')
+
+        for d in SM_HIGG_DECAYS:
+            self.f_BR_scalings.append("CVktkl_BRscal_%s"%d)
+
+    def BuildHscalings(self):
+        # 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")
+
+        for production in SM_HIGG_PROD:
+
+            if production in [ "ggZH", "tHq", "tHW"]: 
+                self.modelBuilder.factory_('expr::CVktkl_pos_XSscal_%s_%s("0.+@0*(@0>0)",Scaling_%s_%s)'\
+                                           %(production,energy,production,energy))                
+                self.f_H_scalings.append("CVktkl_pos_XSscal_%s_%s"%(production,energy))
+
+            elif production in [ "ggH", "qqH" ]:
+                EWK = EWKmap_13[production]
+                self.modelBuilder.factory_('expr::CVktkl_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.modelBuilder.factory_('expr::CVktkl_pos_XSscal_%s_%s("0.+@0*(@0>0)",CVktkl_XSscal_%s_%s)'\
+                                           %(production,energy,production,energy))
+                self.f_H_scalings.append("CVktkl_pos_XSscal_%s_%s"%(production,energy))
+
+            elif production in [ "ZH", "WH", "VH"]: 
+                EWK = EWKmap_13[production]
+                self.modelBuilder.factory_('expr::CVktkl_XSscal_%s_%s("(@1*@1+(@0-1)*%g/%g)/((1-(@0*@0-1)*%g))",kl,CV)'\
+                                           %(production,energy,cXSmap_13[production],EWK,dZH))
+                self.modelBuilder.factory_('expr::CVktkl_pos_XSscal_%s_%s("0.+@0*(@0>0)",CVktkl_XSscal_%s_%s)'\
+                                           %(production,energy,production,energy))
+                self.f_H_scalings.append("CVktkl_pos_XSscal_%s_%s"%(production,energy))
+
+            elif production == "ttH": 
+                EWK = EWKmap_13[production]
+                self.modelBuilder.factory_('expr::CVktkl_XSscal_%s_%s("(@1*@1+(@0-1)*%g/%g)/((1-(@0*@0-1)*%g))",kl,kt)'\
+                                           %(production,energy,cXSmap_13[production],EWK,dZH))
+                self.modelBuilder.factory_('expr::CVktkl_pos_XSscal_%s_%s("0.+@0*(@0>0)",CVktkl_XSscal_%s_%s)'\
+                                           %(production,energy,production,energy))
+                self.f_H_scalings.append("CVktkl_pos_XSscal_%s_%s"%(production,energy))
+
+
+    def findBRscalings(self,process):
+        # process must have the syntax production_whatever_BR1(BR2) 
+        # e.g. ggHH_kl_1_kt_1_2016_hbbhgg, or ttH_2016_Ilikepizza_hzz
+        BRstr = process.split('_')[-1]
+        proc_f_BR_scalings = ['CVktkl_BRscal_h'+x[:2] for x in BRstr.split('h')[1:] ]#[1:] because first substring is always empty, [:2] because we want only the first 3 characters oH decay, e.g. hzz4l --> hzz  
+        # make sure that at least one BR is parsed, and that they are all supported  
+        if len(proc_f_BR_scalings)==0:
+            raise Exception("error parsing branching ratio(s) from {} for signal process {}".format(
+                BRstr, process))
+
+        for proc_f_BR_scaling in proc_f_BR_scalings:
+            if not (proc_f_BR_scaling in self.f_BR_scalings):
+                raise Exception("braching ratio scaling {} for signal process {} is not supported".format(
+                    proc_f_BR_scaling, process))
+                
+        return proc_f_BR_scalings,BRstr
+                
+    def find_singleHmatch(self,process):
+        # single H processes must have the syntax production_whatever(_BR) 
+        # e.g. ttH_2016
+        production = process.split('_')[0]
+        proc_f_H_scaling = "CVktkl_pos_XSscal_%s_%s"%(production,energy)
+        if not (proc_f_H_scaling in self.f_H_scalings):
+            #I am dealing either with an unsupported process or a background process --> I will possibly print a warning
+            if not process in self.nowarning_processes:
+                print "[WARNING] process {} might be not supported. If this is a single or double H process consider investigating.".format(
+                    process)
+                self.nowarning_processes.append(process)#avoid to print multiple times the same message
+            return ""
+        return proc_f_H_scaling
+
+    def BuildXSBRscalingHH(self,f_XS_scaling, process):
+        BRscalings,BRstr = self.findBRscalings(process)
+        if len(BRscalings)!=2:
+            raise Exception("Found {} braching ratio(s) instead of two for signal process {}".format(
+                len(BRscalings), process))
+        f_XSBR_scaling = "%s_BRscal_%s"%(f_XS_scaling,BRstr)
+        if not self.modelBuilder.out.function(f_XSBR_scaling):
+            self.modelBuilder.factory_('expr::%s("@0*@1*@2",%s,%s,%s)'%(f_XSBR_scaling,f_XS_scaling,BRscalings[0],BRscalings[1]))
+        return f_XSBR_scaling
+
+    def BuildXSBRscalingH(self,f_XS_scaling, process):
+        BRscalings,BRstr = self.findBRscalings(process)
+        if len(BRscalings)!=1:
+            raise Exception("Found {} braching ratio(s) instead of one for signal process {}".format(
+                len(BRscalings), process))
+        f_XSBR_scaling = "%s_BRscal_%s"%(f_XS_scaling,BRstr)
+        if not self.modelBuilder.out.function(f_XSBR_scaling):
+            self.modelBuilder.factory_('expr::%s("@0*@1",%s,%s)'%(f_XSBR_scaling,f_XS_scaling,BRscalings[0]))
+
+        return f_XSBR_scaling
+
+