Skip to content

Commit

Permalink
implemented dark matter flux and events calculation for COHERENT expe…
Browse files Browse the repository at this point in the history
…riment.
  • Loading branch information
Ikaroshu committed Jun 10, 2019
1 parent 60c0854 commit e963906
Show file tree
Hide file tree
Showing 2 changed files with 162 additions and 54 deletions.
135 changes: 104 additions & 31 deletions pyCEvNS/events.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,35 +2,12 @@
CEvNS events
"""

from abc import ABC, abstractmethod

from scipy.special import spherical_jn

from .detectors import *
from .flux import *


class EventGen(ABC):
"""
abstract class for events generator,
the inherited class must implement rates and events function
"""
@abstractmethod
def rates(self, er, **kwargs):
"""
:param er: recoil energy in MeV
:return: dN/dE_r
"""
pass

@abstractmethod
def events(self, ea, eb, **kwargs):
"""
:return: number of events in [ea, eb]
"""
pass


def formfsquared(q, rn=5.5, **kwargs):
"""
form factor squared
Expand Down Expand Up @@ -239,7 +216,7 @@ def rates(er):
expo * mev_per_kg * 24 * 60 * 60 / np.dot(det.m, det.frac)


class NSIEventsGen(EventGen):
class NSIEventsGen:
def __init__(self, flux: Flux, detector: Detector, expo: float, target='nucleus', nsi_param=NSIparameters(),
osci_param=oscillation_parameters(), osci_func=None, formfactsq=formfsquared, q2form=False, efficiency=None):
self.flux = flux
Expand Down Expand Up @@ -274,15 +251,111 @@ def events(self, ea, eb, flavor='e', **kwargs):
return Exception('Target should be either nucleus or electron!')


def rates_dm(er, det: Detector, fx: DMFlux, **kwargs):
return np.dot(det.frac, e_charge**4 * fx.epsi_dm**2 * fx.epsi_quark**2 * det.z**2 *
def rates_dm(er, det: Detector, fx: DMFlux, mediator_mass=None, epsilon=None, efficiency=None, **kwargs):
"""
calculating dark matter scattering rates per nucleus
:param er: recoil energy in MeV
:param det: detector
:param fx: dark matter flux
:param mediator_mass: mediator mass in MeV
:param epsilon: mediator to quark coupling multiply by mediator to dark matter coupling
:param efficiency: efficiency function
:return: dark matter scattering rates per nucleus
"""
if mediator_mass is None:
mediator_mass = fx.dp_mass
if epsilon is None:
epsilon = fx.epsi_quark
if efficiency is not None:
return np.dot(det.frac, e_charge**4 * epsilon**2 * det.z**2 *
(2*det.m*fx.fint2(er, det.m) -
(er+(det.m**2*er-fx.dm_mass**2*er)/(2*det.m))*2*det.m*fx.fint(er, det.m) +
er**2*det.m*fx.fint(er, det.m)) / (4*np.pi*(2*det.m*er+fx.dp_mass**2)**2) *
formfsquared(np.sqrt(2*er*det.m), **kwargs))
er**2*det.m*fx.fint(er, det.m)) / (4*np.pi*(2*det.m*er+mediator_mass**2)**2) *
formfsquared(np.sqrt(2*er*det.m), **kwargs)) * efficiency(er)
else:
return np.dot(det.frac, e_charge ** 4 * epsilon ** 2 * det.z ** 2 *
(2 * det.m * fx.fint2(er, det.m) -
(er + (det.m ** 2 * er - fx.dm_mass ** 2 * er) / (2 * det.m)) * 2 * det.m * fx.fint(er, det.m) +
er**2 * det.m * fx.fint(er, det.m)) / (4 * np.pi * (2 * det.m * er + mediator_mass**2)**2) *
formfsquared(np.sqrt(2 * er * det.m), **kwargs))


def binned_events_dm(era, erb, expo, det: Detector, fx: DMFlux, **kwargs):
def binned_events_dm(era, erb, expo, det: Detector, fx: DMFlux, mediator_mass=None, epsilon=None, efficiency=None, **kwargs):
"""
:return: number of nucleus recoil events in the bin [era, erb]
"""
def rates(er):
return rates_dm(er, det, fx, **kwargs)
return quad(rates, era, erb)[0] * expo * mev_per_kg * 24 * 60 * 60 / np.dot(det.m, det.frac)
return rates_dm(er, det, fx, mediator_mass, epsilon, efficiency, **kwargs)
return quad(rates, era, erb)[0] * expo * mev_per_kg * 24 * 60 * 60 / np.dot(det.m, det.frac)


class DmEventsGen:
"""
Dark matter events generator for COHERENT
"""
def __init__(self, dark_photon_mass, life_time, dark_matter_mass, detector_type='csi',
detector_distance=19.3, pot_mu=0.75, pot_sigma=0.25, size=100000):
self.dp_mass = dark_photon_mass
self.tau = life_time
self.dm_mass = dark_matter_mass
self.det_dist = detector_distance
self.mu = pot_mu
self.sigma = pot_sigma
self.size = size
self.det = Detector(detector_type)
self.fx = None
self.generate_flux()

def generate_flux(self):
self.fx = DMFlux(self.dp_mass, self.tau, 1, self.dm_mass, self.det_dist, self.mu, self.sigma, self.size)

def set_dark_photon_mass(self, dark_photon_mass):
self.dp_mass = dark_photon_mass
self.generate_flux()

def set_life_time(self, life_time):
self.tau = life_time
self.generate_flux()

def set_dark_matter_mass(self, dark_matter_mass):
self.dm_mass = dark_matter_mass
self.generate_flux()

def set_detector_distance(self, detector_distance):
self.det_dist = detector_distance
self.generate_flux()

def set_pot_mu(self, pot_mu):
self.mu = pot_mu
self.generate_flux()

def set_pot_sigma(self, pot_sigma):
self.sigma = pot_sigma
self.generate_flux()

def set_size(self, size):
self.size = size
self.generate_flux()

def events(self, mediator_mass, epsilon, n_meas):
"""
generate events according to the time and energy in measured data
:param mediator_mass: mediator mass
:param epsilon: mediator coupling to quark multiply by mediator coupling to dark matter
:param n_meas: measured data
:return: predicted number of event according to the time and energy in the measured data
"""
pe_per_mev = 0.0878 * 13.348 * 1000
n_dm = np.zeros(n_meas.shape[0])
tmin = n_meas[:, 1].min()
plist = np.zeros(int((n_meas[:, 1].max()-tmin)/0.5)+1)
for tt in self.fx.timing:
if int((tt-tmin+0.25)/0.5) < plist.shape[0]:
plist[int((tt-tmin+0.25)/0.5)] += 1
plist /= self.fx.timing.shape[0]
for i in range(n_meas.shape[0]):
pe = n_meas[i, 0]
t = n_meas[i, 1]
n_dm[i] = binned_events_dm((pe - 1)/pe_per_mev, (pe + 1)/pe_per_mev, 4466,
self.det, self.fx, mediator_mass, epsilon, eff_coherent, rn=5.5) * plist[int((t-tmin)/0.5)]
return n_dm
81 changes: 58 additions & 23 deletions pyCEvNS/flux.py
Original file line number Diff line number Diff line change
Expand Up @@ -154,7 +154,7 @@ def fx(ev):
return self.flux(ev, flavor, f, **kwargs)

if not isinstance(emin, np.ndarray):
res = quad(fx, emin, self.evMax)[0]
res = quad(fx, emin, self.evMax)[0] # no need to check range, because outside evMin and evMax are 0
if self.fl_name == 'solar':
if f is None:
f = survival_solar
Expand Down Expand Up @@ -316,6 +316,17 @@ class DMFlux:
"""
def __init__(self, dark_photon_mass, life_time, coupling_quark, dark_matter_mass,
detector_distance=19.3, pot_mu=0.75, pot_sigma=0.25, size=100000):
"""
initialize and generate flux
:param dark_photon_mass: dark photon mass
:param life_time: life time of dark photon in rest frame, unit in micro second
:param coupling_quark: dark photon coupling to quarks
:param dark_matter_mass: mass of dark matter, unit in MeV
:param detector_distance: distance from the detector to the Hg target
:param pot_mu: mean of guassian distribution of proton on target, unit in micro second
:param pot_sigma: std of guassian distribution of proton on target, unit in micro second
:param size: size of sampling dark photons
"""
self.dp_mass = dark_photon_mass
self.dm_mass = dark_matter_mass
self.epsi_quark = coupling_quark
Expand All @@ -324,10 +335,17 @@ def __init__(self, dark_photon_mass, life_time, coupling_quark, dark_matter_mass
self.mu = pot_mu * 1e-6 * c_light / meter_by_mev
self.sigma = pot_sigma * 1e-6 * c_light / meter_by_mev
self.timing, self.energy = self._generate(size)
self.ed_min = self.energy.min()
self.ed_max = self.energy.max()
self.dm_norm = self.epsi_quark**2*0.23*1e20 / (4*np.pi*(detector_distance**2)*24*3600) * (meter_by_mev**2) * \
self.timing.shape[0] * 2 / size

def _generate(self, size=1000000):
"""
generate dark matter flux at COHERENT
:param size: size of sampling dark photons
:return: time and energy histogram of dark matter
"""
dp_m = self.dp_mass
dp_e = ((massofpi+massofp)**2 - massofn**2 + dp_m**2)/(2*(massofpi+massofp))
dp_p = np.sqrt(dp_e ** 2 - dp_m ** 2)
Expand Down Expand Up @@ -364,56 +382,73 @@ def _generate(self, size=1000000):
return np.array(timing) / c_light * meter_by_mev * 1e6, np.array(energy)

def flux(self, ev):
pass
"""
dark matter flux
:param ev: dark matter energy
:return: dark matter flux
"""
return 1/(self.ed_max-self.ed_min)*self.dm_norm if self.ed_min <= ev <= self.ed_max else 0

def fint(self, er, m, **kwargs):
"""
flux/(ex^2-mx^2) integration
:param er:
:param m:
:param kwargs:
:return:
:param er: recoil energy in MeV
:param m: target nucleus mass in MeV
:param kwargs: other argument
:return: flux/(ex^2-mx^2) integration
"""
emin = 0.5 * (np.sqrt((er**2*m+2*er*m**2+2*er*self.dm_mass**2+4*m*self.dm_mass**2)/m) + er)

def integrand(ex):
return self.flux(ex)/(ex**2 - self.dm_mass**2)

if not isinstance(emin, np.ndarray):
res = 1/(self.ex**2-self.dm_mass**2) if emin < self.ex else 0
res = quad(integrand, emin, self.ed_max)[0]
else:
res = np.zeros_like(emin)
for i in range(emin.shape[0]):
res[i] = 1/(self.ex**2-self.dm_mass**2) if emin[i] < self.ex else 0
return res*self.dm_norm
res[i] = quad(integrand, emin[i], self.ed_max)[0]
return res

def fint1(self, er, m, **kwargs):
"""
flux*ex/(ex^2-mx^2) integration
:param er:
:param m:
:param kwargs:
:return:
:param er: recoil energy in MeV
:param m: target nucleus mass in MeV
:param kwargs: other argument
:return: flux*ex/(ex^2-mx^2) integration
"""
emin = 0.5 * (np.sqrt((er**2*m+2*er*m**2+2*er*self.dm_mass**2+4*m*self.dm_mass**2)/m) + er)

def integrand(ex):
return self.flux(ex) * ex / (ex ** 2 - self.dm_mass ** 2)

if not isinstance(emin, np.ndarray):
res = self.ex/(self.ex**2-self.dm_mass**2) if emin < self.ex else 0
res = quad(integrand, emin, self.ed_max)[0]
else:
res = np.zeros_like(emin)
for i in range(emin.shape[0]):
res[i] = self.ex/(self.ex**2-self.dm_mass**2) if emin[i] < self.ex else 0
return res*self.dm_norm
res[i] = quad(integrand, emin[i], self.ed_max)[0]
return res

def fint2(self, er, m, **kwargs):
"""
flux*ex^2/(ex^2-mx^2) integration
:param er:
:param m:
:param kwargs:
:return:
:param er: recoil energy in MeV
:param m: target nucleus mass in MeV
:param kwargs: other argument
:return: flux*ex^2/(ex^2-mx^2) integration
"""
emin = 0.5 * (np.sqrt((er**2*m+2*er*m**2+2*er*self.dm_mass**2+4*m*self.dm_mass**2)/m) + er)

def integrand(ex):
return self.flux(ex) * ex**2 / (ex ** 2 - self.dm_mass ** 2)

if not isinstance(emin, np.ndarray):
res = self.ex**2/(self.ex**2-self.dm_mass**2) if emin < self.ex else 0
res = quad(integrand, emin, self.ed_max)[0]
else:
res = np.zeros_like(emin)
for i in range(emin.shape[0]):
res[i] = self.ex**2/(self.ex**2-self.dm_mass**2) if emin[i] < self.ex else 0
return res*self.dm_norm
res[i] = quad(integrand, emin[i], self.ed_max)[0]
return res

0 comments on commit e963906

Please sign in to comment.