HiggsJPC.py

from HiggsAnalysis.CombinedLimit.PhysicsModel import *


### This base class implements signal yields by production and decay mode
### Specific models can be obtained redefining getHiggsSignalYieldScale
class TwoHypotesisHiggs(PhysicsModel):
    def __init__(self):
        self.mHRange = []
        self.muAsPOI = False
        self.muFloating = False
        self.fqqIncluded = False
        self.fqqFloating = False
        self.poiMap = []
        self.pois = {}
        self.verbose = False
        self.altSignal = "ALT"

    def setModelBuilder(self, modelBuilder):
        PhysicsModel.setModelBuilder(self, modelBuilder)
        self.modelBuilder.doModelBOnly = False

    def getYieldScale(self, bin, process):
        "Split in production and decay, and call getHiggsSignalYieldScale; return 1 for backgrounds"
        if not self.DC.isSignal[process]:
            return 1

        isAlt = self.altSignal in process

        if self.pois:
            target = "%(bin)s/%(process)s" % locals()
            scale = 1
            for p, l in self.poiMap:
                for el in l:
                    if re.match(el, target):
                        scale = p + self.sigNorms[isAlt]

            print("Will scale ", target, " by ", scale)
            return scale

        elif self.fqqIncluded:
            ret = self.sigNorms[isAlt]
            if isAlt:
                ret += self.sigNormsqqH["qqbarH" in process]
            print("Process ", process, " will get scaled by ", ret)
            return ret

        else:
            print("Process ", process, " will get norm ", self.sigNorms[isAlt])
            return self.sigNorms[isAlt]

    def setPhysicsOptions(self, physOptions):
        for po in physOptions:
            if po == "fqqIncluded":
                print("Will consider fqq = fraction of qqH in Alt signal (signal strength will be left floating)")
                # Here alsways setting muFloating if fqq in model, should this be kept optional?
                self.fqqIncluded = True
                self.muFloating = True
            if po == "fqqFloating":
                self.fqqIncluded = True
                self.fqqFloating = True
                self.fqqRange = "0.", "1."
                self.muFloating = True
            if po.startswith("fqqRange="):
                self.fqqIncluded = True
                self.fqqFloating = True
                self.muFloating = True
                self.fqqRange = po.replace("fqqRange=", "").split(",")
                if len(self.fqqRange) != 2:
                    raise RuntimeError("fqq range definition requires two extrema")
                elif float(self.fqqRange[0]) >= float(self.fqqRange[1]):
                    raise RuntimeError("Extrema for fqq range defined with inverterd order. Second must be larger than the first")
            if po == "muAsPOI":
                print("Will consider the signal strength as a parameter of interest")
                self.muAsPOI = True
                self.muFloating = True
            if po == "muFloating":
                print(
                    "Will consider the signal strength as a floating parameter (as a parameter of interest if --PO muAsPOI is specified, as a nuisance otherwise)"
                )
                self.muFloating = True
            if po.startswith("altSignal="):
                self.altSignal = po.split(",")[1]
            if po.startswith("higgsMassRange="):
                self.mHRange = po.replace("higgsMassRange=", "").split(",")
                if len(self.mHRange) != 2:
                    raise RuntimeError("Higgs mass range definition requires two extrema")
                elif float(self.mHRange[0]) >= float(self.mHRange[1]):
                    raise RuntimeError("Extrema for Higgs mass range defined with inverterd order. Second must be larger than the first")
            if po.startswith("verbose"):
                self.verbose = True
            if po.startswith("map="):
                self.muFloating = True
                maplist, poi = po.replace("map=", "").split(":")
                maps = maplist.split(",")
                poiname = re.sub(r"\[.*", "", poi)
                if poiname not in self.pois:
                    if self.verbose:
                        print("Will create a var ", poiname, " with factory ", poi)
                    self.pois[poiname] = poi
                if self.verbose:
                    print("Mapping ", poiname, " to ", maps, " patterns")
                self.poiMap.append((poiname, maps))

    def doParametersOfInterest(self):
        """Create POI and other parameters, and define the POI set."""
        self.modelBuilder.doVar("x[0,0,1]")
        poi = "x"

        if self.muFloating:
            if self.pois:
                for pn, pf in self.pois.items():
                    self.modelBuilder.doVar(pf)
                    if self.muAsPOI:
                        print("Treating %(pn)s as a POI" % locals())
                        poi += "," + pn

                    self.modelBuilder.factory_('expr::%(pn)s_times_not_x("@0*(1-@1)", %(pn)s, x)' % locals())
                    self.modelBuilder.factory_('expr::%(pn)s_times_x("@0*@1", %(pn)s, x)' % locals())
                self.sigNorms = {True: "_times_x", False: "_times_not_x"}

            else:
                self.modelBuilder.doVar("r[1,0,4]")
                if self.muAsPOI:
                    print("Treating r as a POI")
                    poi += ",r"

                self.modelBuilder.factory_('expr::r_times_not_x("@0*(1-@1)", r, x)')
                self.modelBuilder.factory_('expr::r_times_x("@0*@1", r, x)')
                self.sigNorms = {True: "r_times_x", False: "r_times_not_x"}

                if self.fqqIncluded:
                    if self.fqqFloating:
                        self.modelBuilder.doVar(f"fqq[0,{self.fqqRange[0]},{self.fqqRange[1]}]")
                    else:
                        self.modelBuilder.doVar("fqq[0]")
                    self.modelBuilder.factory_('expr::r_times_x_times_fqq("@0*@1", r_times_x, fqq)')
                    self.modelBuilder.factory_('expr::r_times_x_times_not_fqq("@0*(1-@1)", r_times_x, fqq)')

                    self.sigNormsqqH = {True: "_times_fqq", False: "_times_not_fqq"}

        else:
            self.modelBuilder.factory_('expr::not_x("(1-@0)", x)')
            self.sigNorms = {True: "x", False: "not_x"}
        if self.modelBuilder.out.var("MH"):
            var = self.modelBuilder.out.var("MH")
            if len(self.mHRange):
                print(
                    "MH will be left floating within",
                    self.mHRange[0],
                    "and",
                    self.mHRange[1],
                )
                var.setRange(float(self.mHRange[0]), float(self.mHRange[1]))
                var.setConstant(False)
                poi += ",MH"
            else:
                print("MH will be assumed to be", self.options.mass)
                var.removeMin()
                var.removeMax()
                var.setVal(self.options.mass)
        else:
            if len(self.mHRange):
                print(
                    "MH will be left floating within",
                    self.mHRange[0],
                    "and",
                    self.mHRange[1],
                )
                self.modelBuilder.doVar(f"MH[{self.mHRange[0]},{self.mHRange[1]}]")
                poi += ",MH"
            else:
                print("MH (not there before) will be assumed to be", self.options.mass)
                self.modelBuilder.doVar("MH[%g]" % self.options.mass)
        self.modelBuilder.doSet("POI", poi)


twoHypothesisHiggs = TwoHypotesisHiggs()