From bba39cc7601a03c51f5d2d4ede350675d689047c Mon Sep 17 00:00:00 2001 From: Shu Date: Thu, 18 Jul 2019 11:15:20 -0500 Subject: [PATCH] fixed survival probability for earth, update a dark matter events generator for COHERENT experiment, can choose color for the plot --- pyCEvNS/events.py | 51 ++++++++++++++++++++++++++-------------- pyCEvNS/oscillation.py | 53 ++++++++++++++++++++++++++++++++++++++---- pyCEvNS/plot.py | 4 ++-- 3 files changed, 84 insertions(+), 24 deletions(-) diff --git a/pyCEvNS/events.py b/pyCEvNS/events.py index ce70f78..914eb52 100644 --- a/pyCEvNS/events.py +++ b/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,26 +267,37 @@ 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) @@ -293,8 +305,8 @@ 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 diff --git a/pyCEvNS/oscillation.py b/pyCEvNS/oscillation.py index e8b5851..a0249af 100644 --- a/pyCEvNS/oscillation.py +++ b/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)) diff --git a/pyCEvNS/plot.py b/pyCEvNS/plot.py index 8143e55..4d1fe1b 100644 --- a/pyCEvNS/plot.py +++ b/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]]), "*")