plotBR.py

import math
import numpy as np
import ROOT as rt

mu = 2.3*0.001
md = 4.8*0.001
mc = 1.275
ms = 95*0.001
mt = 173.2
mb = 4.18
me = 0.000511
mmu = 0.1066
mtau = 1.777
Nc=3.
alphas=0.1177
v=246.

def gamma_scalar_bb(mphi, mchi, gq, gl, gchi):

    aux0=(mb**2)*(mphi*(Nc*((np.maximum(0.,(1.+(-4.*((mb**2)*(mphi**-2.))))))**1.5)))
    output=((1./484128.)*((gq**2)*aux0))/math.pi
    
    return output

def gamma_scalar(mphi, mchi, gq, gl, gchi):

    # quarks
    aux0=(mu**2)*(mphi*(Nc*((np.maximum(0.,(1.+(-4.*((mu**2)*(mphi**-2.))))))**1.5)))
    aux1=(md**2)*(mphi*(Nc*((np.maximum(0.,(1.+(-4.*((md**2)*(mphi**-2.))))))**1.5)))
    aux2=(mc**2)*(mphi*(Nc*((np.maximum(0.,(1.+(-4.*((mc**2)*(mphi**-2.))))))**1.5)))
    aux3=(ms**2)*(mphi*(Nc*((np.maximum(0.,(1.+(-4.*((ms**2)*(mphi**-2.))))))**1.5)))
    aux4=(mt**2)*(mphi*(Nc*((np.maximum(0.,(1.+(-4.*((mt**2)*(mphi**-2.))))))**1.5)))
    aux5=(mb**2)*(mphi*(Nc*((np.maximum(0.,(1.+(-4.*((mb**2)*(mphi**-2.))))))**1.5)))

    output0=((1./484128.)*((gq**2)*aux0))/math.pi
    output1=((1./484128.)*((gq**2)*aux1))/math.pi
    output2=((1./484128.)*((gq**2)*aux2))/math.pi
    output3=((1./484128.)*((gq**2)*aux3))/math.pi
    output4=((1./484128.)*((gq**2)*aux4))/math.pi
    output5=((1./484128.)*((gq**2)*aux5))/math.pi

    # leptons
    aux6=(me**2)*(mphi*((np.maximum(0.,(1.+(-4.*((me**2)*(mphi**-2.))))))**1.5))
    aux7=(mmu**2)*(mphi*((np.maximum(0.,(1.+(-4.*((mmu**2)*(mphi**-2.))))))**1.5))
    aux8=(mtau**2)*(mphi*((np.maximum(0.,(1.+(-4.*((mtau**2)*(mphi**-2.))))))**1.5))
    
    output6=((1./484128.)*((gl**2)*aux6))/math.pi
    output7=((1./484128.)*((gl**2)*aux7))/math.pi
    output8=((1./484128.)*((gl**2)*aux8))/math.pi

    # dark matter
    aux9=(gchi**2)*(mphi*((np.maximum(0.,(1.+(-4.*((mchi**2)*(mphi**-2.))))))**1.5))
    output9=(0.125*aux9)/math.pi
    
    # gluons (loop)
    aux10=(1.+(-4.*((mphi**-2.)*(mt**2))))*(((np.arctan((1./np.sqrt(-1.+(4.*((mphi**-2.)*(mt**2)))+0j))))**2))
    aux11=(math.pi**-3.)*((v**-2.)*(((np.absolute(((mphi**-2.)*((mt**2)*(1.+aux10)))))**2)));
    output10=0.0000165246*((alphas**2)*((gq**2)*((mphi**3.)*((mt**2)*aux11))));
    
    return output0+output1+output2+output3+output4+output5+output6+output7+output8+output9+output10

def br_scalar_bb_mphi(x, par):
    mphi = x[0]
    mchi = par[0]
    gq = par[1]
    gl = par[2]
    gchi = par[3]

    return gamma_scalar_bb(mphi, mchi, gq, gl, gchi)/gamma_scalar(mphi, mchi, gq, gl, gchi)

def gamma_pseudoscalar_bb(mphi, mchi, gq, gl, gchi):

    aux0=(mb**2)*(mphi*(Nc*(np.sqrt((np.maximum(0.,(1.+(-4.*((mb**2)*(mphi**-2.))))))))))
    output=((1./484128.)*((gq**2)*aux0))/math.pi
    
    return output

def gamma_pseudoscalar(mphi, mchi, gq, gl, gchi):

    # quarks
    aux0=(mu**2)*(mphi*(Nc*(np.sqrt((np.maximum(0.,(1.+(-4.*((mu**2)*(mphi**-2.))))))))))
    aux1=(md**2)*(mphi*(Nc*(np.sqrt((np.maximum(0.,(1.+(-4.*((md**2)*(mphi**-2.))))))))))
    aux2=(mc**2)*(mphi*(Nc*(np.sqrt((np.maximum(0.,(1.+(-4.*((mc**2)*(mphi**-2.))))))))))
    aux3=(ms**2)*(mphi*(Nc*(np.sqrt((np.maximum(0.,(1.+(-4.*((ms**2)*(mphi**-2.))))))))))
    aux4=(mt**2)*(mphi*(Nc*(np.sqrt((np.maximum(0.,(1.+(-4.*((mt**2)*(mphi**-2.))))))))))
    aux5=(mb**2)*(mphi*(Nc*(np.sqrt((np.maximum(0.,(1.+(-4.*((mb**2)*(mphi**-2.))))))))))

    output0=((1./484128.)*((gq**2)*aux0))/math.pi
    output1=((1./484128.)*((gq**2)*aux1))/math.pi
    output2=((1./484128.)*((gq**2)*aux2))/math.pi
    output3=((1./484128.)*((gq**2)*aux3))/math.pi
    output4=((1./484128.)*((gq**2)*aux4))/math.pi
    output5=((1./484128.)*((gq**2)*aux5))/math.pi

    # leptons
    aux6=(me**2)*(mphi*(np.sqrt((np.maximum(0.,(1.+(-4.*((me**2)*(mphi**-2.)))))))))
    aux7=(mmu**2)*(mphi*(np.sqrt((np.maximum(0.,(1.+(-4.*((mmu**2)*(mphi**-2.)))))))))
    aux8=(mtau**2)*(mphi*(np.sqrt((np.maximum(0.,(1.+(-4.*((mtau**2)*(mphi**-2.)))))))))
    
    output6=((1./484128.)*((gl**2)*aux6))/math.pi
    output7=((1./484128.)*((gl**2)*aux7))/math.pi
    output8=((1./484128.)*((gl**2)*aux8))/math.pi

    # dark matter    
    aux9=(gchi**2)*(mphi*(np.sqrt((np.maximum(0.,(1.+(-4.*((mchi**2)*(mphi**-2.)))))))))
    output9=(0.125*aux9)/math.pi
    
    # gluons (loop)
    if mphi == 2*mt:
        aux10=math.pi**2/4.
    else:
        aux10=(mphi**-2.)*((mt**2)*(((np.arctan((1./np.sqrt(-1.+(4.*((mphi**-2.)*(mt**2)))+0j))))**2)))
        
    aux11=(gq**2)*((mphi**3.)*((mt**2)*((math.pi**-3.)*((v**-2.)*(((np.absolute(aux10))**2))))))
    output10=0.0000165246*((alphas**2)*aux11)
    
    return output0+output1+output2+output3+output4+output5+output6+output7+output8+output9+output10

def br_pseudoscalar_bb_mphi(x, par):
    mphi = x[0]
    mchi = par[0]
    gq = par[1]
    gl = par[2]
    gchi = par[3]

    return gamma_pseudoscalar_bb(mphi, mchi, gq, gl, gchi)/gamma_pseudoscalar(mphi, mchi, gq, gl, gchi)


rt.gStyle.SetOptStat(0)
rt.gStyle.SetOptTitle(0)
c= rt.TCanvas('c','c',500,400)
hist = rt.TH1D('hist','hist',100,0,600)
hist.SetMaximum(1)
hist.SetMinimum(0)
hist.GetXaxis().SetRangeUser(0, 600)

f1 = rt.TF1('br_DMSbb',br_pseudoscalar_bb_mphi, 0, 600, 4)
f1.SetParameter(0, 1.) #mchi
f1.SetParameter(1, 1.) #gq
f1.SetParameter(2, 0.) #gl
f1.SetParameter(3, 1.) #gchi


f2 = rt.TF1('br_DMPSbb',br_pseudoscalar_bb_mphi, 0, 600, 4)
f2.SetParameter(0, 1.) #mchi
f2.SetParameter(1, 1.) #gq
f2.SetParameter(2, 0.) #gl
f2.SetParameter(3, 1.) #gchi

hist.Draw()
f1.Draw("csame")
c.Print('br_scalar_bb_mphi.pdf')

hist.Draw()
f2.Draw("csame")
c.Print('br_pseudoscalar_bb_mphi.pdf')