proton_bremsstrahlung.py

Go to the documentation of this file.

# SPDX-License-Identifier: LGPL-3.0-or-later
# SPDX-FileCopyrightText: Copyright CERN for the benefit of the SHiP Collaboration
 
import math
import sys
 
import numpy as np
import ROOT as r
from darkphoton import alphaQED
from scipy.integrate import dblquad
 
# proton mass
mProton = 0.938272081  # GeV/c - PDG2016
protonEnergy = 400.0  # GeV/c
protonMomentum = math.sqrt(protonEnergy * protonEnergy - mProton * mProton)
 
 
# VDM FORM FACTOR
def rhoFormFactor(m):
    """From https://arxiv.org/abs/0910.5589"""
    # constants from the code from Inar: https://github.com/pgdeniverville/BdNMC/blob/master/src/Proton_Brem_Distribution.cpp
    f1ra = 0.6165340033101271
    f1rb = 0.22320420111672623
    f1rc = -0.33973820442685326
    f1wa = 1.0117544786579074
    f1wb = -0.8816565944110686
    f1wc = 0.3699021157531611
    f1prho = f1ra * 0.77**2 / (0.77**2 - m**2 - 0.77 * 0.15j)
    f1prhop = f1rb * 1.25**2 / (1.25**2 - m**2 - 1.25 * 0.3j)
    f1prhopp = f1rc * 1.45**2 / (1.45**2 - m**2 - 1.45 * 0.5j)
    f1pomega = f1wa * 0.77**2 / (0.77**2 - m**2 - 0.77 * 0.0085j)
    f1pomegap = f1wb * 1.25**2 / (1.25**2 - m**2 - 1.25 * 0.3j)
    f1pomegapp = f1wc * 1.45**2 / (1.45**2 - m**2 - 1.45 * 0.5j)
    return abs(f1prho + f1prhop + f1prhopp + f1pomega + f1pomegap + f1pomegapp)
 
 
# useful functions
def energy(p, m):
    """Compute energy from momentum and mass"""
    return math.sqrt(p * p + m * m)
 
 
def penaltyFactor(m):
    """Penalty factor for high masses - dipole form factor in the proton-A' vertex"""
    """ m in GeV """
    if m * m > 0.71:
        return math.pow(m * m / 0.71, -4)
    else:
        return 1
 
 
def zeta(p, theta):
    """Fraction of the proton momentum carried away by the paraphoton in the beam direction"""
    return p / (protonMomentum * math.sqrt(theta * theta + 1.0))
 
 
def pTransverse(p, theta):
    """Paraphoton transverse momentum in the lab frame"""
    return protonMomentum * theta * zeta(p, theta)
 
 
def ptSquare(p, theta):
    """Square paraphoton transverse momentum in the lab frame"""
    return pow(pTransverse(p, theta), 2.0)
 
 
def H(p, theta, mDarkPhoton):
    """A kinematic term"""
    return (
        ptSquare(p, theta)
        + (1.0 - zeta(p, theta)) * mDarkPhoton * mDarkPhoton
        + pow(zeta(p, theta), 2.0) * mProton * mProton
    )
 
 
def wba(p, theta, mDarkPhoton, epsilon):
    """Cross section weighting function in the Fermi-Weizsaeker-Williams approximation"""
    const = epsilon * epsilon * alphaQED / (2.0 * math.pi * H(p, theta, mDarkPhoton))
 
    h2 = pow(H(p, theta, mDarkPhoton), 2.0)
    oneMinusZSquare = pow(1.0 - zeta(p, theta), 2.0)
    mp2 = mProton * mProton
    mA2 = mDarkPhoton * mDarkPhoton
 
    p1 = (1.0 + oneMinusZSquare) / zeta(p, theta)
    p2 = (
        2.0
        * zeta(p, theta)
        * (1.0 - zeta(p, theta))
        * ((2.0 * mp2 + mA2) / H(p, theta, mDarkPhoton) - pow(zeta(p, theta), 2.0) * 2.0 * mp2 * mp2 / h2)
    )
    # p3 = 2.*zeta(p,theta)*(1.-zeta(p,theta))*(zeta(p,theta)+oneMinusZSquare)*mp2*mA2/h2
    p3 = 2.0 * zeta(p, theta) * (1.0 - zeta(p, theta)) * (1 + oneMinusZSquare) * mp2 * mA2 / h2
    p4 = 2.0 * zeta(p, theta) * oneMinusZSquare * mA2 * mA2 / h2
    return const * (p1 - p2 + p3 + p4)
 
 
def sigma(s):  # s in GeV^2 ---> sigma in mb
    """Parametrisation of sigma(s)"""
    a1 = 35.45
    a2 = 0.308
    a3 = 28.94
    a4 = 33.34
    a5 = 0.545
    a6 = 0.458
    a7 = 42.53
    p1 = a2 * pow(math.log(s / a3), 2.0)
    p2 = a4 * pow((1.0 / s), a5)
    p3 = a7 * pow((1.0 / s), a6)
    return a1 + p1 - p2 + p3
 
 
def es(p, mDarkPhoton):
    """s(p,mA)"""
    return 2.0 * mProton * (energy(protonMomentum, mProton) - energy(p, mDarkPhoton))
 
 
def sigmaRatio(p, mDarkPhoton):
    """sigma(s') / sigma(s)"""
    return sigma(es(p, mDarkPhoton)) / sigma(2.0 * mProton * energy(protonMomentum, mProton))
 
 
def dNdZdPtSquare(p, mDarkPhoton, theta, epsilon):
    """Differential A' rate per p.o.t. as a function of Z and Pt^2"""
    return sigmaRatio(p, mDarkPhoton) * wba(p, theta, mDarkPhoton, epsilon)
 
 
def dPt2dTheta(p, theta):
    """Jacobian Pt^2->theta"""
    z2 = pow(zeta(p, theta), 2.0)
    return 2.0 * theta * z2 * protonMomentum * protonMomentum
 
 
def dZdP(p, theta):
    """Jacobian z->p"""
    return 1.0 / (protonMomentum * math.sqrt(theta * theta + 1.0))
 
 
def dNdPdTheta(p, theta, mDarkPhoton, epsilon):
    """Differential A' rate per p.o.t. as a function of P and theta"""
    diffRate = dNdZdPtSquare(p, mDarkPhoton, theta, epsilon) * dPt2dTheta(p, theta) * dZdP(p, theta)
    return math.fabs(diffRate)  # integrating in (-pi, pi)...
 
 
def pMin(mDarkPhoton):
    return max(0.14 * protonMomentum, mDarkPhoton)
 
 
def pMax(mDarkPhoton):
    # return min(0.86*protonMomentum, math.sqrt( (energy(protonMomentum,mProton)**2. - mDarkPhoton**2.) - mDarkPhoton**2.))
    return math.sqrt((energy(protonMomentum, mProton) ** 2.0 - mDarkPhoton**2.0) - mDarkPhoton**2.0)
 
 
def prodRate(mDarkPhoton, epsilon, tmin=-0.5 * math.pi, tmax=0.5 * math.pi):
    """dNdPdTheta integrated over p and theta"""
    integral = dblquad(
        dNdPdTheta,  # integrand
        tmin,
        tmax,  # theta boundaries (2nd argument of integrand)
        lambda x: pMin(mDarkPhoton),
        lambda x: pMax(mDarkPhoton),  # p boundaries (1st argument of integrand)
        args=(mDarkPhoton, epsilon),
    )  # extra parameters to pass to integrand
    return integral[0]
 
 
# total production rate of A'
# norm = prodRate(1.1,3.e-7) #mDarkPhoton,epsilon)
# number of A' produced
# numDarkPhotons = int(math.floor(norm*protonFlux))
#
# print
# print "Epsilon \t %s"%epsilon
# print "mDarkPhoton \t %s"%mDarkPhoton
# print "A' production rate per p.o.t: \t %.8g"%norm
# print "Number of A' produced in SHiP: \t %.8g"%numDarkPhotons
 
 
def normalisedProductionPDF(p, theta, mDarkPhoton, epsilon, norm):
    """Probability density function for A' production in SHIP"""
    return (1.0 / norm) * dNdPdTheta(p, theta, mDarkPhoton, epsilon)
 
 
def hProdPDF(mDarkPhoton, epsilon, norm, binsp, binstheta, tmin=-0.5 * math.pi, tmax=0.5 * math.pi, suffix=""):
    """Histogram of the PDF for A' production in SHIP"""
    angles = np.linspace(tmin, tmax, binstheta).tolist()
    anglestep = 2.0 * (tmax - tmin) / binstheta
    momentumStep = (pMax(mDarkPhoton) - pMin(mDarkPhoton)) / (binsp - 1)
    momenta = np.linspace(pMin(mDarkPhoton), pMax(mDarkPhoton), binsp, endpoint=False).tolist()
    hPDF = r.TH2F(
        f"hPDF_eps{epsilon}_m{mDarkPhoton}",
        f"hPDF_eps{epsilon}_m{mDarkPhoton}",
        binsp,
        pMin(mDarkPhoton) - 0.5 * momentumStep,
        pMax(mDarkPhoton) - 0.5 * momentumStep,
        binstheta,
        tmin - 0.5 * anglestep,
        tmax - 0.5 * anglestep,
    )
    hPDF.SetTitle(f"PDF for A' production (m_{{A'}}={mDarkPhoton} GeV, #epsilon ={epsilon})")
    hPDF.GetXaxis().SetTitle("P_{A'} [GeV]")
    hPDF.GetYaxis().SetTitle("#theta_{A'} [rad]")
    hPDFtheta = r.TH1F(
        f"hPDFtheta_eps{epsilon}_m{mDarkPhoton}",
        f"hPDFtheta_eps{epsilon}_m{mDarkPhoton}",
        binstheta,
        tmin - 0.5 * anglestep,
        tmax - 0.5 * anglestep,
    )
    hPDFp = r.TH1F(
        f"hPDFp_eps{epsilon}_m{mDarkPhoton}",
        f"hPDFp_eps{epsilon}_m{mDarkPhoton}",
        binsp,
        pMin(mDarkPhoton) - 0.5 * momentumStep,
        pMax(mDarkPhoton) - 0.5 * momentumStep,
    )
    hPDFp.GetXaxis().SetTitle("P_{A'} [GeV]")
    hPDFtheta.GetXaxis().SetTitle("#theta_{A'} [rad]")
    for theta in angles:
        for p in momenta:
            w = normalisedProductionPDF(p, theta, mDarkPhoton, epsilon, norm)
            hPDF.Fill(p, theta, w)
            hPDFtheta.Fill(theta, w)
            hPDFp.Fill(p, w)
    hPdfFilename = sys.modules["__main__"].options.outputDir + f"/ParaPhoton_eps{epsilon}_m{mDarkPhoton}{suffix}.root"
    outfile = r.TFile(hPdfFilename, "recreate")
    # weight = hPDF.Integral("width")
    # print "Weight = %3.3f"%weight
    # hPDF.Scale(1./weight)
    hPDF.Write()
    hPDFp.Write()
    hPDFtheta.Write()
    outfile.Close()
    del angles
    del momenta
    return hPDF