Reproducing the COHERENT results - and New Physics constraints

Code for reproducing the CEvNS signal observed by COHERENT - see arXiv:1708.01294. Note that the COHERENT-2017 data are now publicly available (arXiv:1804.09459) - this notebook uses digitized results from the original 2017 paper.

Note that we neglect the axial charge of the nucleus, and thus the contribution from strange quarks. We also use a slightly different parametrisation of the Form Factor, compared to the COHERENT collaboration.



In [1]:

    
from __future__ import print_function

%matplotlib inline
import numpy as np
import matplotlib
#matplotlib.use('Agg')
import matplotlib.pyplot as pl

from scipy.integrate import quad
from scipy.interpolate import interp1d, UnivariateSpline,InterpolatedUnivariateSpline

from scipy.optimize import minimize

from tqdm import tqdm

#Change default font size so you don't need a magnifying glass
matplotlib.rc('font', **{'size'   : 16})

Import the CEvNS module (for calculating the signal spectrum and loading the neutrino fluxes)



In [2]:

    
import CEvNS
#help(CEvNS.xsec_CEvNS)

Neutrino Flux @ SNS

Let's load the neutrino flux. Note that here we're only plotting the continuum. There is also a population of monochromatic (29.65 MeV) muon neutrinos which we add in separately in the code (because the flux is a delta-function, it's hard to model here).



In [3]:

    
#Initialise neutrino_flux interpolation function
CEvNS.loadNeutrinoFlux("SNS")

#Plot neutrino flux
E_nu = np.logspace(0, np.log10(300),1000)

pl.figure()
pl.semilogy(E_nu, CEvNS.neutrino_flux_tot(E_nu))

pl.title(r"Neutrino flux at SNS", fontsize=12)
pl.xlabel(r"Neutrino energy, $E_\nu$ [MeV]")
pl.ylabel(r"$\Phi_\nu$ [cm$^{-2}$ s$^{-1}$ MeV$^{-1}$]")
pl.show()









    



/Users/bradkav/anaconda3/lib/python3.6/site-packages/scipy/integrate/quadpack.py:364: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  warnings.warn(msg, IntegrationWarning)

COHERENT efficiency function

Load in the efficiency (as a function of photoelectrons, PE). Set to zero below 5 PE.



In [4]:

    
COHERENT_PE, COHERENT_eff = np.loadtxt("DataFiles/COHERENT_eff.txt", unpack=True)
effinterp = interp1d(COHERENT_PE, COHERENT_eff, bounds_error=False, fill_value=0.0)

def efficiency_single(x):
    if (x > 4.9):
        return effinterp(x)
    else:
        return 1e-10
    
efficiency = np.vectorize(efficiency_single)

PEvals = np.linspace(0, 50, 100)
pl.figure()
pl.plot(PEvals, efficiency(PEvals))
pl.xlabel("PE")
pl.ylabel("Efficiency")
pl.show()

COHERENT event rate

Calculate number of CEvNS signal events at COHERENT (in bins of 2 PE)



In [5]:

    
#Nuclear properties for Cs and I
A_Cs = 133.0
Z_Cs = 55.0

A_I = 127.0
Z_I = 53.0

#Mass fractions
f_Cs = A_Cs/(A_Cs + A_I)
f_I = A_I/(A_Cs + A_I)

mass = 14.6 #target mass in kg
time = 308.1 #exposure time in days

PEperkeV = 1.17 #Number of PE per keV

#Get the differential rate function from the CEvNS module
#Note that this function allows for an extra vector mediator, 
#but the default coupling is zero, so we'll forget about it
diffRate_CEvNS = CEvNS.differentialRate_CEvNS
    
#Differential rates (times efficiency) for the two target nuclei, per PE
dRdPE_Cs = lambda x: (1.0/PEperkeV)*efficiency(x)*mass*time*f_Cs*diffRate_CEvNS(x/PEperkeV, A_Cs, Z_Cs)
dRdPE_I = lambda x: (1.0/PEperkeV)*efficiency(x)*mass*time*f_I*diffRate_CEvNS(x/PEperkeV, A_I, Z_I)



In [6]:

    
#Calculate number of signal events in each bin in the Standard Model (SM)
PE_bins = np.linspace(0, 50, 26)

N_SM_Cs = np.zeros(25)
N_SM_I = np.zeros(25)
N_SM_tot = np.zeros(25)

for i in tqdm(range(25)):
    N_SM_Cs[i] = quad(dRdPE_Cs, PE_bins[i], PE_bins[i+1], epsabs = 0.01)[0]
    N_SM_I[i] = quad(dRdPE_I, PE_bins[i], PE_bins[i+1], epsabs = 0.01)[0]
    N_SM_tot[i] = N_SM_Cs[i] + N_SM_I[i]
    
print("Total CEvNS events expected: ", np.sum(N_SM_tot))









    



100%|██████████| 25/25 [00:06<00:00,  3.82it/s]





    



Total CEvNS events expected:  167.179814234625

Comparing with the COHERENT results

First, let's load in the observed data and calculated spectrum (digitized from arXiv:1708.01294).



In [7]:

    
COHERENT_data = np.loadtxt("DataFiles/COHERENT_data.txt", usecols=(1,))
COHERENT_upper = np.loadtxt("DataFiles/COHERENT_upper.txt", usecols=(1,)) - COHERENT_data
COHERENT_lower = COHERENT_data - np.loadtxt("DataFiles/COHERENT_lower.txt", usecols=(1,))
COHERENT_spect = np.loadtxt("DataFiles/COHERENT_spectrum.txt", usecols=(1,))
COHERENT_bins = np.arange(1,50,2)

Now plot the results:



In [8]:

    
pl.figure(figsize=(10,6))

pl.step(PE_bins, np.append(N_SM_tot,0), 'g', linestyle="-", where = "post", label="CEvNS signal (this work)",linewidth=1.5)
pl.step(PE_bins, np.append(COHERENT_spect,0), 'g', linestyle="--", where = "post", label="CEvNS signal (1708.01294)",linewidth=1.5)

pl.axhline(0, linestyle='--', color = 'gray')

pl.errorbar(COHERENT_bins, COHERENT_data, fmt='ko', \
            yerr = [COHERENT_lower, COHERENT_upper], label="COHERENT data",\
            capsize=0.0)

pl.xlabel("Number of photoelectrons (PE)")
pl.ylabel("Res. counts / 2 PE")

pl.legend( fontsize=14)

pl.xlim(0, 50)
pl.ylim(-15, 35)

pl.savefig("plots/COHERENT_data.pdf", bbox_inches="tight")
pl.show()

Fit to signal strength

Very simple fit to the number of CEvNS signal events, using only a 1-bin likelihood.

We start by defining the $\chi^2$, as given in arXiv:1708.01294. We use a generic form, so that we don't have to recalculate the number of signal events all the time...



In [9]:

    
def chisq_generic(N_sig, alpha, beta):   

    
    #Beam-on backgrounds
    N_BG = 6.0
    #Number of measured events
    N_meas = 142.0 
    #Statistical uncertainty
    sig_stat = np.sqrt(N_meas + 2*405 + N_BG) 
    
    #Uncertainties
    unc = (alpha/0.28)**2 + (beta/0.25)**2
    
    return ((N_meas - N_sig*(1.0+alpha) - N_BG*(1.0+beta))**2)/sig_stat**2 + unc



In [10]:

    
#Calculate minimum chi-squared as a function of (alpha, beta) nuisance parameters
def minchisq_Nsig(Nsig):
    
    minres = minimize(lambda x: chisq_generic(Nsig, x[0], x[1]), (0.0,0.0))
    return minres.fun



In [11]:

    
Nsiglist= np.linspace(0, 1000,1001)
chi2list = [minchisq_Nsig(Ns) for Ns in Nsiglist]

delta_chi2 = (chi2list - np.min(chi2list))



In [12]:

    
pl.figure(figsize=(6,6))
pl.plot(Nsiglist, delta_chi2, linewidth=2.0)
pl.ylim(0, 25)

pl.axvline(np.sum(N_SM_tot), linestyle='--', color='k')
pl.text(172, 20, "SM prediction")

pl.ylabel(r"$\Delta \chi^2$")
pl.xlabel(r"CE$\nu$NS counts")

pl.savefig("plots/COHERENT_likelihood.pdf", bbox_inches="tight")

pl.show()

To speed things up later (so we don't have to do the minimization every time), we'll tabulate and interpolate the chi-squared as a function of the number of signal events. This works because we're using a simple chi-squared which depends only on the number of signal events:



In [13]:

    
deltachi2_Nsig = interp1d(Nsiglist, delta_chi2, bounds_error=False, fill_value=delta_chi2[-1])

NSI constraints

Calculate constraints on NSI parameters. Here, we're just assuming that the flavor-conserving e-e NSI couplings are non-zero, so we have to calculate the contribution to the rate from only the electron neutrinos and then see how that changes:



In [14]:

    
#Differential rates (times efficiency) for the two target nuclei, per PE 
# For electron neutrinos ONLY
dRdPE_Cs_e = lambda x: (1.0/PEperkeV)*efficiency(x)*mass*time*f_Cs*diffRate_CEvNS(x/PEperkeV, A_Cs, Z_Cs, nu_flavor="e")
dRdPE_I_e = lambda x: (1.0/PEperkeV)*efficiency(x)*mass*time*f_I*diffRate_CEvNS(x/PEperkeV, A_I, Z_I, nu_flavor="e")

# For muon neutrinos ONLY
dRdPE_Cs_mu = lambda x: (1.0/PEperkeV)*efficiency(x)*mass*time*f_Cs*(diffRate_CEvNS(x/PEperkeV, A_Cs, Z_Cs, nu_flavor="mu")+ diffRate_CEvNS(x/PEperkeV, A_Cs, Z_Cs, nu_flavor="mub"))
dRdPE_I_mu = lambda x: (1.0/PEperkeV)*efficiency(x)*mass*time*f_I*(diffRate_CEvNS(x/PEperkeV, A_I, Z_I, nu_flavor="mu") + diffRate_CEvNS(x/PEperkeV, A_I, Z_I, nu_flavor="mub"))


#Now calculate bin-by-bin signal from electron neutrinos
bins_Cs_e = np.zeros(25)
bins_I_e = np.zeros(25)
for i in tqdm(range(25)):
    bins_Cs_e[i] = quad(dRdPE_Cs_e, PE_bins[i], PE_bins[i+1], epsabs = 0.01)[0]
    bins_I_e[i] = quad(dRdPE_I_e, PE_bins[i], PE_bins[i+1], epsabs = 0.01)[0]

print("Number of CEvNS events due to nu_e: ", np.sum(bins_Cs_e + bins_I_e))

#Now calculate bin-by-bin signal from muon neutrinos
bins_Cs_mu = np.zeros(25)
bins_I_mu = np.zeros(25)
for i in tqdm(range(25)):
    bins_Cs_mu[i] = quad(dRdPE_Cs_mu, PE_bins[i], PE_bins[i+1], epsabs = 0.01)[0]
    bins_I_mu[i] = quad(dRdPE_I_mu, PE_bins[i], PE_bins[i+1], epsabs = 0.01)[0]

print("Number of CEvNS events due to nu_mu: ", np.sum(bins_Cs_mu + bins_I_mu))









    



100%|██████████| 25/25 [00:06<00:00,  3.97it/s]
  4%|▍         | 1/25 [00:00<00:03,  7.67it/s]





    



Number of CEvNS events due to nu_e:  86.3440584916171






    



100%|██████████| 25/25 [00:12<00:00,  1.93it/s]





    



Number of CEvNS events due to nu_mu:  80.87450071444458

Flavour-conserving NSI

Now, let's calculate the correction to the CEvNS rate from flavor-conserving NSI:



In [15]:

    
def NSI_corr(eps_uV, eps_dV, A, Z):
    SIN2THETAW = 0.2387
    
    #Calculate standard weak nuclear charge (squared)
    Qsq = 4.0*((A - Z)*(-0.5) + Z*(0.5 - 2*SIN2THETAW))**2
    
    #Calculate the modified nuclear charge from NSI
    Qsq_NSI = 4.0*((A - Z)*(-0.5 + eps_uV + 2.0*eps_dV) + Z*(0.5 - 2*SIN2THETAW + 2*eps_uV + eps_dV))**2
    
    return Qsq_NSI/Qsq

Calculate simplified (single bin) chi-squared (see chi-squared expression around p.32 in COHERENT paper):



In [16]:

    
def deltachisq_NSI_ee(eps_uV, eps_dV):
    #NB: bins_I and bins_Cs are calculated further up in the script (they are the SM signal prediction)

    #Signal events from Iodine (with NSI correction only applying to electron neutrino events)
    N_sig_I = (N_SM_I + (NSI_corr(eps_uV, eps_dV, A_I, Z_I) - 1.0)*bins_I_e)
    #Now signal events from Caesium
    N_sig_Cs = (N_SM_Cs + (NSI_corr(eps_uV, eps_dV, A_Cs, Z_Cs) - 1.0)*bins_Cs_e)
    
    #Number of signal events
    N_NSI = np.sum(N_sig_I + N_sig_Cs)
    
    return deltachi2_Nsig(N_NSI)

Calculate the (minimum) chi-squared on a grid and save to file:



In [17]:

    
Ngrid = 101

ulist = np.linspace(-1.0, 1.0, Ngrid)
dlist = np.linspace(-1.0, 1.0, Ngrid)
UL, DL = np.meshgrid(ulist, dlist)
delta_chi2_grid_ee = 0.0*UL


#Not very elegant loop
for i in tqdm(range(Ngrid)):
    for j in range(Ngrid):
        delta_chi2_grid_ee[i,j] = deltachisq_NSI_ee(UL[i,j], DL[i,j])

#Find best-fit point
ind_BF = np.argmin(delta_chi2_grid_ee)
BF = [UL.flatten()[ind_BF], DL.flatten()[ind_BF]]
        
print("Best fit point: ", BF)
np.savetxt("results/COHERENT_NSI_deltachi2_ee.txt", delta_chi2_grid_ee, header="101x101 grid, corresponding to (uV, dV) values between -1 and 1. Flavor-conserving ee NSI.")









    



100%|██████████| 101/101 [00:00<00:00, 226.29it/s]





    



Best fit point:  [-0.76, 1.0]

Plot the 90% allowed regions:



In [18]:

    
pl.figure(figsize=(6,6))
#pl.contourf(DL, UL, delta_chi2_grid, levels=[0,1,2,3,4,5,6,7,8,9,10],cmap="Blues")
pl.contourf(DL, UL, delta_chi2_grid_ee, levels=[0,4.6],cmap="Blues")
#levels=[0,4.60]
#pl.colorbar()

pl.plot(0.0, 0.0,'k+', markersize=12.0, label="Standard Model")

pl.plot(BF[1], BF[0], 'ro', label="Best fit")
#pl.plot(-0.25, 0.5, 'ro')

pl.ylabel(r"$\epsilon_{ee}^{uV}$", fontsize=22.0)
pl.xlabel(r"$\epsilon_{ee}^{dV}$" ,fontsize=22.0)

pl.title(r"$90\%$ CL allowed regions", fontsize=16.0)

pl.legend(frameon=False, fontsize=12, numpoints=1)

pl.savefig("plots/COHERENT_NSI_ee.pdf", bbox_inches="tight")
pl.show()

Flavour-changing NSI ($e\mu$)

Now the correction to the CEvNS rate from flavor-changing NSI ($e\mu$-type):



In [19]:

    
def NSI_corr_changing(eps_uV, eps_dV, A, Z):
    SIN2THETAW = 0.2387
    
    #Calculate standard weak nuclear charge (squared)
    Qsq = 4.0*((A - Z)*(-0.5) + Z*(0.5 - 2*SIN2THETAW))**2
    
    #Calculate the modified nuclear charge from NSI
    Qsq_NSI = Qsq + 4.0*((A-Z)*(eps_uV + 2.0*eps_dV) + Z*(2.0*eps_uV + eps_dV))**2
    
    return Qsq_NSI/Qsq



In [20]:

    
def deltachisq_NSI_emu(eps_uV, eps_dV):
    #NB: bins_I and bins_Cs are calculated further up in the script (they are the SM signal prediction)

    N_sig_I = (N_SM_I)*NSI_corr_changing(eps_uV, eps_dV, A_I, Z_I)
    #Now signal events from Caesium
    N_sig_Cs = (N_SM_Cs)*NSI_corr_changing(eps_uV, eps_dV, A_Cs, Z_Cs)
    
    #Number of signal events
    N_NSI = np.sum(N_sig_I + N_sig_Cs)
    
    return deltachi2_Nsig(N_NSI)

Calculate delta-chisquared over a grid and save to file



In [21]:

    
Ngrid = 101

ulist = np.linspace(-1.0, 1.0, Ngrid)
dlist = np.linspace(-1.0, 1.0, Ngrid)
UL, DL = np.meshgrid(ulist, dlist)
delta_chi2_grid_emu = 0.0*UL

#Not very elegant loop
for i in tqdm(range(Ngrid)):
    for j in range(Ngrid):
        delta_chi2_grid_emu[i,j] = deltachisq_NSI_emu(UL[i,j], DL[i,j])

#Find best-fit point
ind_BF = np.argmin(delta_chi2_grid_emu)
BF = [UL.flatten()[ind_BF], DL.flatten()[ind_BF]]
        
print("Best fit point: ", BF)
np.savetxt("results/COHERENT_NSI_deltachi2_emu.txt", delta_chi2_grid_emu, header="101x101 grid, corresponding to (uV, dV) values between -1 and 1.")









    



100%|██████████| 101/101 [00:00<00:00, 262.78it/s]





    



Best fit point:  [0.0, 0.0]



In [22]:

    
pl.figure(figsize=(6,6))
#pl.contourf(DL, UL, delta_chi2_grid, levels=[0,1,2,3,4,5,6,7,8,9,10],cmap="Blues")
pl.contourf(DL, UL, delta_chi2_grid_emu, levels=[0,4.6],cmap="Blues")
#levels=[0,4.60]
#pl.colorbar()

pl.plot(0.0, 0.0,'k+', markersize=12.0, label="Standard Model")

pl.plot(BF[1], BF[0], 'ro', label="Best fit")
#pl.plot(-0.25, 0.5, 'ro')

pl.ylabel(r"$\epsilon_{e\mu}^{uV}$", fontsize=22.0)
pl.xlabel(r"$\epsilon_{e\mu}^{dV}$" ,fontsize=22.0)

pl.title(r"$90\%$ CL allowed regions", fontsize=16.0)

pl.legend(frameon=False, fontsize=12, numpoints=1)

pl.savefig("plots/COHERENT_NSI_emu.pdf", bbox_inches="tight")
pl.show()

Flavour-changing NSI ($e\tau$)

Finally, allowed regions for Flavour-changing NSI ($e\tau$-type)



In [23]:

    
def deltachisq_NSI_etau(eps_uV, eps_dV):
    #NB: bins_I and bins_Cs are calculated further up in the script (they are the SM signal prediction)

    #Signal events from Iodine (with NSI correction only applying to electron neutrino events)
    N_sig_I = (N_SM_I + (NSI_corr_changing(eps_uV, eps_dV, A_I, Z_I) - 1.0)*bins_I_e)
    #Now signal events from Caesium
    N_sig_Cs = (N_SM_Cs + (NSI_corr_changing(eps_uV, eps_dV, A_Cs, Z_Cs) - 1.0)*bins_Cs_e)
    
    #Number of signal events
    N_NSI = np.sum(N_sig_I + N_sig_Cs)
    
    return deltachi2_Nsig(N_NSI)



In [24]:

    
Ngrid = 101

ulist = np.linspace(-1.0, 1.0, Ngrid)
dlist = np.linspace(-1.0, 1.0, Ngrid)
UL, DL = np.meshgrid(ulist, dlist)
delta_chi2_grid_etau = 0.0*UL

#Not very elegant loop
for i in tqdm(range(Ngrid)):
    for j in range(Ngrid):
        delta_chi2_grid_etau[i,j] = deltachisq_NSI_etau(UL[i,j], DL[i,j])

#Find best-fit point
ind_BF = np.argmin(delta_chi2_grid_etau)
BF = [UL.flatten()[ind_BF], DL.flatten()[ind_BF]]
        
print("Best fit point: ", BF)
np.savetxt("results/COHERENT_NSI_deltachi2_etau.txt", delta_chi2_grid_etau, header="101x101 grid, corresponding to (uV, dV) values between -1 and 1.")









    



100%|██████████| 101/101 [00:00<00:00, 233.46it/s]





    



Best fit point:  [0.0, 0.0]



In [25]:

    
pl.figure(figsize=(6,6))
#pl.contourf(DL, UL, delta_chi2_grid, levels=[0,1,2,3,4,5,6,7,8,9,10],cmap="Blues")
pl.contourf(DL, UL, delta_chi2_grid_etau, levels=[0,4.6],cmap="Blues")
#levels=[0,4.60]
#pl.colorbar()

pl.plot(0.0, 0.0,'k+', markersize=12.0, label="Standard Model")

pl.plot(BF[1], BF[0], 'ro', label="Best fit")
#pl.plot(-0.25, 0.5, 'ro')

pl.ylabel(r"$\epsilon_{e\tau}^{uV}$", fontsize=22.0)
pl.xlabel(r"$\epsilon_{e\tau}^{dV}$" ,fontsize=22.0)

pl.title(r"$90\%$ CL allowed regions", fontsize=16.0)

pl.legend(frameon=False, fontsize=12, numpoints=1)

pl.savefig("plots/COHERENT_NSI_etau.pdf", bbox_inches="tight")
pl.show()

Limits on the neutrino magnetic moment

Now let's calculate a limit on the neutrino magnetic moment (again, from a crude single-bin $\chi^2$).



In [26]:

    
#Calculate the number of neutrino magnetic moment scattering events
#assuming a universal magnetic moment (in units of 1e-12 mu_B)

diffRate_mag = np.vectorize(CEvNS.differentialRate_magnetic)

dRdPE_mag = lambda x: (1.0/PEperkeV)*efficiency(x)*mass*time*(f_Cs*diffRate_mag(x/PEperkeV, A_Cs, Z_Cs, 1e-12)\
                                                              + f_I*diffRate_mag(x/PEperkeV, A_I, Z_I, 1e-12))
N_mag = quad(dRdPE_mag, 0, 50)[0]
print("Number of magnetic moment signal events (for mu_nu = 1e-12 mu_B):", N_mag)









    



Number of magnetic moment signal events (for mu_nu = 1e-12 mu_B): 9.239043774475988e-06



In [27]:

    
def deltachisq_mag(mu_nu):    
    
    #Signal events is sum of standard CEvNS + magnetic moment events
    N_sig = np.sum(N_SM_tot) + N_mag*(mu_nu/1e-12)**2
    
    return deltachi2_Nsig(N_sig)

Scan over a grid:



In [28]:

    
Ngrid = 501

maglist = np.logspace(-12, -6, Ngrid)



deltachi2_list_mag = 0.0*maglist
#Not very elegant loop
for i in tqdm(range(Ngrid)):
    deltachi2_list_mag[i] = deltachisq_mag(maglist[i])

upper_limit = maglist[deltachi2_list_mag > 2.706][0]

print("90% upper limit: ", upper_limit)









    



100%|██████████| 501/501 [00:00<00:00, 27277.46it/s]





    



90% upper limit:  3.3728730865886923e-09

Do some plotting:



In [29]:

    
pl.figure(figsize=(6,6))
pl.semilogx(maglist, deltachi2_list_mag, linewidth=2.0)
#pl.ylim(0, 25)

pl.axhline(2.706, linestyle='--', color='k')
pl.axvline(upper_limit, linestyle=':', color='k')

pl.text(1e-11, 3, "90% CL")

pl.ylabel(r"$\Delta \chi^2$")
pl.xlabel(r"Neutrino magnetic moment, $\mu_{\nu} / \mu_B$")

pl.savefig("plots/COHERENT_magnetic.pdf", bbox_inches="tight")

pl.show()

Limits on new vector mediators

First, let's calculate the total number of signal events at a given mediator mass and coupling...

It takes a while to recalculate the number of signal events for each mediator mass and coupling, so we'll do some rescaling and interpolation trickery:



In [30]:

    
def tabulate_rate( m_med):
    vector_rate = lambda x, gsq: (1.0/PEperkeV)*efficiency(x)*mass*time*(f_Cs*CEvNS.differentialRate_CEvNS(x/PEperkeV, A_Cs, Z_Cs,gsq,m_med)\
                                                              + f_I*CEvNS.differentialRate_CEvNS(x/PEperkeV, A_I, Z_I, gsq,m_med))
    alpha = 1.0
    
    PE_min = 4.0
    PE_max = 50.0
    Nvals = 500
    PEvals = np.logspace(np.log10(PE_min), np.log10(PE_max),Nvals)
    
    Rvals_A = [np.sqrt(vector_rate(PEvals[i], 0)) for i in range(Nvals)]
    Rvals_B = [(1.0/(4.0*alpha*Rvals_A[i]))*(vector_rate(PEvals[i], alpha) - vector_rate(PEvals[i], -alpha)) for i in range(Nvals)]
    
    tabrate_A = InterpolatedUnivariateSpline(PEvals, Rvals_A, k = 1)
    tabrate_B = InterpolatedUnivariateSpline(PEvals, Rvals_B, k = 1)
    
    return tabrate_A, tabrate_B



In [31]:

    
def N_sig_vector(gsq, m_med):
    integrand = lambda x: (1.0/PEperkeV)*efficiency(x)*mass*time*(f_Cs*CEvNS.differentialRate_CEvNS(x/PEperkeV, A_Cs, Z_Cs,gsq,m_med)\
                                                              + f_I*CEvNS.differentialRate_CEvNS(x/PEperkeV, A_I, Z_I, gsq,m_med))

    xlist = np.linspace(4,50,100)
    integ_vals =  np.vectorize(integrand)(xlist)
    
    return np.trapz(integ_vals, xlist)

def N_sig_vector_tab(gsq, tabrate_A, tabrate_B):
    integrand = lambda x: (tabrate_A(x) + tabrate_B(x)*gsq)**2.0
    
    xlist = np.linspace(4,50,100)
    integ_vals =  np.vectorize(integrand)(xlist)
    
    return np.trapz(integ_vals, xlist)
    
    #return quad(integrand, 4.0, 50, epsabs=0.01)[0]



In [32]:

    
def tabulate_Nsig(tabrate_A, tabrate_B):
    N_A = N_sig_vector_tab(0, tabrate_A, tabrate_B)
    N_C = 0.5*(N_sig_vector_tab(1.0, tabrate_A, tabrate_B) + N_sig_vector_tab(-1.0, tabrate_A, tabrate_B))- N_A
    N_B = N_sig_vector_tab(1.0, tabrate_A, tabrate_B) - N_A - N_C
    return N_A, N_B, N_C



In [33]:

    
def N_sig_fulltab(gsq, Nsig_A, Nsig_B, Nsig_C):
    return Nsig_A + gsq*Nsig_B + gsq**2*Nsig_C



In [34]:

    
#Calculate the number of signal events for a 1000 MeV Z', with coupling 1e-4 by doing:
rate_A, rate_B = tabulate_rate(1000)
N_A, N_B,N_C = tabulate_Nsig(rate_A, rate_B)
#N_sig_vector_tab(1e-4, rate_A, rate_B)
N_sig_fulltab(1e-4, N_A, N_B, N_C)









    Out[34]:





666862.515674649

Now we scan over a grid in $g^2$ and $m_V$ to calculate the $\chi^2$ at each point:



In [35]:

    
gsq_list = np.append(np.logspace(0, 2, 100),1e20)
m_list = np.sort(np.append(np.logspace(-2, 4,49), [1e-6,1e8]))

#Need to search for the limit in a narrow band of coupling values
g_upper = 1e-11*(50**2+m_list**2)
g_lower = 1e-13*(50**2+m_list**2)

deltachi2_vec_grid = np.zeros((51, 101))

for i in tqdm(range(len(m_list))):
    rate_A, rate_B = tabulate_rate(m_list[i])
    N_A, N_B,N_C = tabulate_Nsig(rate_A, rate_B)
    for j, gsq in enumerate(gsq_list):
        N_sig = N_sig_fulltab(gsq*g_lower[i], N_A, N_B, N_C)
        deltachi2_vec_grid[i, j] = deltachi2_Nsig(N_sig)









    



100%|██████████| 51/51 [04:40<00:00,  5.50s/it]



In [36]:

    
mgrid, ggrid = np.meshgrid(m_list, gsq_list, indexing='ij')
ggrid *= 1e-13*(50**2 + mgrid**2)

np.savetxt("results/COHERENT_Zprime.txt", np.c_[mgrid.flatten(), ggrid.flatten(), deltachi2_vec_grid.flatten()])

pl.figure(figsize=(6,6))
pl.loglog(m_list, g_upper, 'k--')
pl.loglog(m_list, g_lower, 'k--')
pl.contourf(mgrid, ggrid, deltachi2_vec_grid, levels=[2.7,1e10],cmap="Blues")

pl.ylim(1e-10, 1e5)
#pl.colorbar()
pl.xlabel(r"$m_{Z'}$ [MeV]")
pl.ylabel(r"$g_{Z'}^2$")
pl.title("Blue region (and above) is excluded...", fontsize=12)
pl.savefig("plots/COHERENT_Zprime.pdf")
pl.show()

Limits on a new scalar mediator

Finally, let's look at limits on the couplings of a new scalar mediator $\phi$. We start by calculating the contribution to the number of signal events for a given mediator mass (this can be rescaled by the coupling $g_\phi^4$ later):



In [37]:

    
def calc_Nsig_scalar(m_med):
    scalar_rate = lambda x: (1.0/PEperkeV)*efficiency(x)*mass*time*(f_Cs*CEvNS.differentialRate_scalar(x/PEperkeV, A_Cs, Z_Cs,1,m_med)\
                                                              + f_I*CEvNS.differentialRate_scalar(x/PEperkeV, A_I, Z_I, 1,m_med))
    
    xlist = np.linspace(4,50,100)
    integ_vals =  np.vectorize(scalar_rate)(xlist)
    
    return np.trapz(integ_vals, xlist)
    
    #return quad(scalar_rate, PE_min, PE_max)[0]

Now grid-scan to get the $\Delta \chi^2$:



In [38]:

    
m_list = np.logspace(-3, 7,50)
gsq_list = np.logspace(0, 4, 50)

#Again, need to search in a specific range of coupling values to find the limit...
g_upper = 1e-10*(50**2+m_list**2)
g_lower = 1e-14*(50**2+m_list**2)

deltachi2_scal_grid = np.zeros((len(m_list), len(gsq_list)))


for i in tqdm(range(len(m_list))):
    Nsig_scalar = calc_Nsig_scalar(m_list[i])
    for j in range(len(gsq_list)):
        deltachi2_scal_grid[i,j] = deltachi2_Nsig(np.sum(N_SM_tot) + Nsig_scalar*(gsq_list[j]*g_lower[i])**2)









    



100%|██████████| 50/50 [00:16<00:00,  3.11it/s]



In [39]:

    
mgrid, ggrid = np.meshgrid(m_list, gsq_list, indexing='ij')
ggrid *= 1e-14*(50**2+mgrid**2)

np.savetxt("results/COHERENT_scalar.txt", np.c_[mgrid.flatten(), ggrid.flatten(), deltachi2_scal_grid.flatten()])


pl.figure(figsize=(6,6))
pl.loglog(m_list, g_upper, 'k--')
pl.loglog(m_list, g_lower, 'k--')
pl.contourf(mgrid, ggrid, deltachi2_scal_grid, levels=[2.7,1e10],cmap="Blues")
#pl.colorbar()
pl.xlabel(r"$m_{\phi}$ [MeV]")
pl.ylabel(r"$g_{\phi}^2$")
pl.title("Blue region (and above) is excluded...", fontsize=12)
pl.savefig("plots/COHERENT_scalar.pdf")
pl.show()



In [ ]: