fixed survival probability for earth,

update a dark matter events generator for COHERENT experiment,
can choose color for the plot

pyCEvNS/events.py

@@ -6,6 +6,7 @@
 
 from .detectors import *
 from .flux import *
+from .helper import _poisson
 
 
 def formfsquared(q, rn=5.5, **kwargs):
@@ -251,7 +252,7 @@ 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, mediator_mass=None, epsilon=None, efficiency=None,  **kwargs):
+def rates_dm(er, det: Detector, fx: DMFlux, mediator_mass=None, epsilon=None, efficiency=None, smear=False,  **kwargs):
     """
     calculating dark matter scattering rates per nucleus
     :param er: recoil energy in MeV
@@ -266,35 +267,46 @@ def rates_dm(er, det: Detector, fx: DMFlux, mediator_mass=None, epsilon=None, ef
         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+mediator_mass**2)**2) *
-                  formfsquared(np.sqrt(2*er*det.m), **kwargs)) * efficiency(er)
+
+    def rates(err):
+        if efficiency is not None:
+            return np.dot(det.frac, e_charge**4 * epsilon**2 * det.z**2 *
+                      (2*det.m*fx.fint2(err, det.m) -
+                       (err+(det.m**2*err-fx.dm_mass**2*err)/(2*det.m))*2*det.m*fx.fint(err, det.m) +
+                       err**2*det.m*fx.fint(err, det.m)) / (4*np.pi*(2*det.m*err+mediator_mass**2)**2) *
+                      formfsquared(np.sqrt(2*err*det.m), **kwargs)) * efficiency(err)
+        else:
+            return np.dot(det.frac, e_charge ** 4 * epsilon ** 2 * det.z ** 2 *
+                          (2 * det.m * fx.fint2(err, det.m) -
+                           (err + (det.m ** 2 * err - fx.dm_mass ** 2 * err) / (2 * det.m)) * 2 * det.m * fx.fint(err, det.m) +
+                           err**2 * det.m * fx.fint(err, det.m)) / (4 * np.pi * (2 * det.m * err + mediator_mass**2)**2) *
+                          formfsquared(np.sqrt(2 * err * det.m), **kwargs))
+
+    if not smear:
+        return rates(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 func(pep):
+            pe_per_mev = 0.0878 * 13.348 * 1000
+            return rates(pep/pe_per_mev) * _poisson(er*pe_per_mev, pep)
+        return quad(func, 0, 60)[0]
 
 
-def binned_events_dm(era, erb, expo, det: Detector, fx: DMFlux, mediator_mass=None, epsilon=None, efficiency=None, **kwargs):
+
+def binned_events_dm(era, erb, expo, det: Detector, fx: DMFlux, mediator_mass=None, epsilon=None, efficiency=None, smear=False, **kwargs):
     """
     :return: number of nucleus recoil events in the bin [era, erb]
     """
     def rates(er):
-        return rates_dm(er, det, fx, mediator_mass, epsilon, efficiency, **kwargs)
+        return rates_dm(er, det, fx, mediator_mass, epsilon, efficiency, smear, **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):
+    def __init__(self, dark_photon_mass, life_time, dark_matter_mass, expo=4466, detector_type='csi',
+                 detector_distance=19.3, pot_mu=0.75, pot_sigma=0.25, size=100000, smear=False, rn=5.5):
         self.dp_mass = dark_photon_mass
         self.tau = life_time
         self.dm_mass = dark_matter_mass
@@ -304,6 +316,9 @@ def __init__(self, dark_photon_mass, life_time, dark_matter_mass, detector_type=
         self.size = size
         self.det = Detector(detector_type)
         self.fx = None
+        self.expo = expo
+        self.smear = smear
+        self.rn = rn
         self.generate_flux()
 
     def generate_flux(self):
@@ -356,6 +371,6 @@ def events(self, mediator_mass, epsilon, n_meas):
         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)]
+            n_dm[i] = binned_events_dm((pe - 1)/pe_per_mev, (pe + 1)/pe_per_mev, self.expo,
+                                       self.det, self.fx, mediator_mass, epsilon, eff_coherent, self.smear, rn=self.rn) * plist[int((t-tmin)/0.5)]
         return n_dm

pyCEvNS/oscillation.py

@@ -137,6 +137,51 @@ def survival_const(ev, lenth=0.0, epsi=NSIparameters(), op=oscillation_parameter
     return np.real(res)
 
 
+def survival_const_amp(ev, lenth=0.0, epsi=NSIparameters(), op=oscillation_parameters(),
+                   ne=2.2 * 6.02e23 * (100 * meter_by_mev) ** 3, nui='e', nuf='e'):
+    """
+    survival/transitional amplitude with constant matter density
+    :param ev: nuetrino energy in MeV
+    :param lenth: oscillation lenth in meters
+    :param epsi: epsilons
+    :param nui: initail flavor
+    :param nuf: final flavor
+    :param op: oscillation parameters
+    :param ne: electron number density in MeV^3
+    :return: survival/transitional probability
+    """
+    dic = {'e': 0, 'mu': 1, 'tau': 2, 'ebar': 0, 'mubar': 1, 'taubar': 2}
+    fi = dic[nui]
+    ff = dic[nuf]
+    lenth = lenth / meter_by_mev
+    opt = op.copy()
+    if nuf[-1] == 'r':
+        opt['delta'] = -opt['delta']
+    o23 = np.array([[1, 0, 0],
+                   [0, np.cos(opt['t23']), np.sin(opt['t23'])],
+                   [0, -np.sin(opt['t23']), np.cos(opt['t23'])]])
+    u13 = np.array([[np.cos(opt['t13']), 0, np.sin(opt['t13']) * (np.exp(- opt['delta'] * 1j))],
+                   [0, 1, 0],
+                   [-np.sin(opt['t13'] * (np.exp(opt['delta'] * 1j))), 0, np.cos(opt['t13'])]])
+    o12 = np.array([[np.cos(opt['t12']), np.sin(opt['t12']), 0],
+                   [-np.sin(opt['t12']), np.cos(opt['t12']), 0],
+                   [0, 0, 1]])
+    umix = o23 @ u13 @ o12
+    m = np.diag(np.array([0, opt['d21'] / (2 * ev), opt['d31'] / (2 * ev)]))
+    vf = np.sqrt(2) * gf * ne * (epsi.ee() + 3 * epsi.eu() + 3 * epsi.ed())
+    if nuf[-1] == 'r':
+        hf = umix @ m @ umix.conj().T - np.conj(vf)
+    else:
+        hf = umix @ m @ umix.conj().T + vf
+    w, v = np.linalg.eigh(hf)
+    res = 0.0
+    for i in range(3):
+        # for j in range(3):
+        theta = (w[i]) * lenth
+        res += v[ff, i] * np.conj(v[fi, i]) * (np.cos(theta) - 1j * np.sin(theta))
+    return res
+
+
 def survival_average(ev, epsi=NSIparameters(), op=oscillation_parameters(),
                      ne=2.2 * 6.02e23 * (100 * meter_by_mev) ** 3, nui='e', nuf='e'):
     opt = op.copy()
@@ -203,7 +248,7 @@ def survial_atmos(ev, zenith=1.0, epsi=NSIparameters(), op=oscillation_parameter
             f_list = ['e', 'mu', 'tau']
         for i in f_list:
             for j in f_list:
-                res += survival_const(ev, l_mantle_half, epsi=epsi, nui=nui, nuf=i, op=op, ne=n_mantle) * \
-                       survival_const(ev, l_core, epsi=epsi, nui=i, nuf=j, ne=n_core) * \
-                       survival_const(ev, l_mantle_half, epsi=epsi, nui=j, nuf=nuf, ne=n_mantle)
-        return res
+                res += survival_const_amp(ev, l_mantle_half, epsi=epsi, nui=nui, nuf=i, op=op, ne=n_mantle) * \
+                       survival_const_amp(ev, l_core, epsi=epsi, nui=i, nuf=j, ne=n_core) * \
+                       survival_const_amp(ev, l_mantle_half, epsi=epsi, nui=j, nuf=nuf, ne=n_mantle)
+        return np.real(res * np.conj(res))

pyCEvNS/plot.py

@@ -85,7 +85,7 @@ def credible_1d(self, idx: int, credible_level=(0.6827, 0.9545), nbins=80, ax=No
         return fig, ax
 
     def credible_2d(self, idx: tuple, credible_level=(0.6827, 0.9545), nbins=80, ax=None,
-                    center=None, heat=False, xlim=None, ylim=None, mark_best=False):
+                    center=None, heat=False, xlim=None, ylim=None, mark_best=False, color='b'):
         """
         plot the correlation between parameters
         :param idx: the index of the two parameters to be ploted
@@ -150,7 +150,7 @@ def credible_2d(self, idx: tuple, credible_level=(0.6827, 0.9545), nbins=80, ax=
                     if s > cl[ic]:
                         cll = zv[sorted_idx[0][i], sorted_idx[1][i]]
                         break
-                ax.contourf(xv, yv, zv, (cll, 1), colors=('b', 'white'), alpha=al[ic])
+                ax.contourf(xv, yv, zv, (cll, 1), colors=(color, 'white'), alpha=al[ic])
             ax.axis('scaled')
         if center is not None:
             ax.plot(np.array([center[0]]), np.array([center[1]]), "*")