Skip to content

Commit c9c677c

Browse files
athompson-gitathompson-git
committed
updating changes with BEBC limit code, improvements in flux calculation
1 parent 2d815fb commit c9c677c

6 files changed

+474
-45
lines changed

Diff for: ‎hnl_flux.py

+9
Original file line numberDiff line numberDiff line change
@@ -345,6 +345,7 @@ def __init__(self, proton_energy=120000.0, target=Material("C"), det_dist=DUNE_D
345345
self.hnl_energy = []
346346
self.hnl_angle = []
347347
self.hnl_flux = []
348+
self.hnl_timing = []
348349
self.decay_axion_weight = []
349350
self.scatter_axion_weight = []
350351
self.support = np.ones(n_samples)
@@ -370,6 +371,7 @@ def set_new_params(self, zprime_mass=None, coupling_BL=None, mixing_angle=None,
370371
self.hnl_energy = []
371372
self.hnl_angle = []
372373
self.hnl_flux = []
374+
self.hnl_timing = []
373375
self.decay_axion_weight = []
374376
self.scatter_axion_weight = []
375377

@@ -380,6 +382,7 @@ def simulate(self):
380382
self.hnl_energy = []
381383
self.hnl_angle = []
382384
self.hnl_flux = []
385+
self.hnl_timing = []
383386
self.decay_axion_weight = []
384387
self.scatter_axion_weight = []
385388

@@ -451,6 +454,7 @@ def propagate(self, Ualpha, timing_cut=None, is_isotropic=False, verbose=False):
451454
tau = boost / gamma_total if gamma_total > 0.0 else np.inf * np.ones_like(boost)
452455

453456
# Calculate time of flight
457+
self.hnl_timing = self.det_dist / (v_a * 1e-2 * C_LIGHT)
454458
if timing_cut is not None:
455459
tof = self.det_dist / (v_a * 1e-2 * C_LIGHT)
456460
delta_tof = abs((self.det_dist / (1e-2 * C_LIGHT)) - tof)
@@ -506,6 +510,7 @@ def __init__(self, meson_flux=[[0.0, 0.0, 0.0, M_PI0]], flux_weight=1.0, meson_s
506510
self.hnl_energy = []
507511
self.hnl_angle = []
508512
self.hnl_flux = []
513+
self.hnl_timing = []
509514
self.decay_axion_weight = []
510515
self.scatter_axion_weight = []
511516
self.support = np.ones(n_samples)
@@ -537,6 +542,7 @@ def set_new_params(self, zprime_mass=None, coupling_BL=None, mixing_angle=None,
537542
self.hnl_energy = []
538543
self.hnl_angle = []
539544
self.hnl_flux = []
545+
self.hnl_timing = []
540546
self.decay_axion_weight = []
541547
self.scatter_axion_weight = []
542548

@@ -556,6 +562,7 @@ def simulate(self):
556562
self.hnl_energy = []
557563
self.hnl_angle = []
558564
self.hnl_flux = []
565+
self.hnl_timing = []
559566
self.decay_axion_weight = []
560567
self.scatter_axion_weight = []
561568

@@ -623,6 +630,7 @@ def decay_to_hnl(self):
623630

624631
def propagate(self, Ualpha, timing_cut=None, is_isotropic=False, verbose=False):
625632

633+
626634
gamma_vis = decay_width_visible(self.hnl_mass, Ualpha, [self.alpha])
627635
gamma_total = total_decay_width_hnl(self.hnl_mass, Ualpha, self.alpha)
628636

@@ -640,6 +648,7 @@ def propagate(self, Ualpha, timing_cut=None, is_isotropic=False, verbose=False):
640648
tau = boost / gamma_total if gamma_total > 0.0 else np.inf * np.ones_like(boost)
641649

642650
# Calculate time of flight
651+
self.hnl_timing = self.det_dist / (v_a * 1e-2 * C_LIGHT)
643652
if timing_cut is not None:
644653
tof = self.det_dist / (v_a * 1e-2 * C_LIGHT)
645654
delta_tof = abs((self.det_dist / (1e-2 * C_LIGHT)) - tof)

Diff for: ‎other_constants.py

+34
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
import sys
2+
sys.path.append("../")
3+
4+
from alplib.constants import *
5+
from alplib.materials import *
6+
7+
# all distances with respect to BNB
8+
9+
# BEBC
10+
# see e.g. https://arxiv.org/pdf/2208.00416
11+
# volume = 3.57 × 2.52 × 1.85
12+
BEBC_POT = 2.72e18 # not certain
13+
BEBC_DIST = 404.0
14+
BEBC_AREA = 3.57*2.52
15+
BEBC_LENGTH = 1.85
16+
BEBC_THRESH = 1000.0
17+
BEBC_EFF = 0.96
18+
BEBC_ET_CUT = 1000.0
19+
# equivalent solid angle: ~ 18 mrad
20+
21+
# NA62
22+
23+
# CHARM
24+
25+
# T2K to ND280
26+
# see https://arxiv.org/pdf/1902.07598
27+
T2K_POT = 12.34e20 # not certain
28+
T2K_DIST = 280.0
29+
T2K_AREA = 3.4
30+
T2K_LENGTH = 1.84
31+
T2K_THRESH = 1000.0
32+
T2K_EFF = 0.20
33+
T2K_ANGLE = 0.0436
34+

Diff for: ‎plot_limits.py

+48-22
Original file line numberDiff line numberDiff line change
@@ -166,25 +166,30 @@ def plot_muon_limits(mass_ratio=2.1):
166166
c1 = get_Usq_upper_lower_limits("scans/hnl_DUNE_scan_mass-ratio-2_ebrem_umu_highstats.txt", color='royalblue')
167167
c2 = get_Usq_contour_from_file("scans/hnl_DUNE_scan_mass-ratio-2_pbrem_umu.txt", color='indianred')
168168
c3 = get_Usq_contour_from_file("scans/hnl_DUNE_scan_mass-ratio-2_meson_umu.txt", color='g')
169+
170+
#c4 = get_Usq_contour_from_file("scans/recast_scans/hnl_BEBC_scan_mass-ratio-2_pbrem_umu.txt", color='k')
171+
169172
elif mass_ratio == 5:
170-
c3 = get_Usq_contour_from_file("scans/hnl_DUNE_scan_mass-ratio-5_pbrem_umu.txt", color='indianred')
171-
c1 = get_Usq_upper_lower_limits("scans/hnl_DUNE_scan_mass-ratio-5_ebrem_umu.txt", color='royalblue')
173+
c1 = get_Usq_contour_from_file("scans/hnl_DUNE_scan_mass-ratio-5_pbrem_umu.txt", color='indianred')
174+
c2 = get_Usq_upper_lower_limits("scans/hnl_DUNE_scan_mass-ratio-5_ebrem_umu.txt", color='royalblue')
172175
c3 = get_Usq_contour_from_file("scans/hnl_DUNE_scan_mass-ratio-5_meson_umu.txt", color='g')
173176

177+
#c4 = get_Usq_contour_from_file("scans/recast_scans/hnl_BEBC_scan_mass-ratio-5_pbrem_umu.txt", color='k')
178+
174179

175180
# plot seesaw line
176181
masses = np.logspace(0, 4, 100)
177182
plt.fill_between(masses, y1=seesaw_lower(masses), y2=seesaw_upper(masses), color='gold', alpha=0.3)
178183

179184
# future dune
180185
plt.plot(1e3*power(10,dune_future[:,0]), power(10, dune_future[:,1]), ls='dotted', color='teal', linewidth=1.0)
181-
plt.text(43.0, 1.2e-8, "DUNE-ND\n(Mixing Only)", rotation=-35.0, color='teal')
186+
plt.text(63.0, 2e-8, "DUNE-ND\n(Mixing Only)", rotation=-32.0, color='teal')
182187

183-
plt.text(30.0, 1e-9, "Seesaw", rotation=-14.0, fontsize=12)
188+
plt.text(100.0, 5e-10, "Seesaw", rotation=-12.0, fontsize=12)
184189
plt.text(1500.0, 5e-4, "Laboratory Limits\n(Mixing-only)", fontsize=12)
185190

186-
line_pbrem = Line2D([0], [0], label=r'$p N \to p N Z^\prime$', color='indianred')
187-
line_ebrem = Line2D([0], [0], label=r'$e^\pm N \to e^\pm N Z^\prime$, $e^+ e^- \to Z^\prime$', color='royalblue')
191+
line_pbrem = Line2D([0], [0], label=r'$p p \to X Z^\prime$', color='indianred')
192+
line_ebrem = Line2D([0], [0], label=r'$e^\pm A \to e^\pm A Z^\prime$, $e^+ e^- \to Z^\prime$', color='royalblue')
188193
line_pion = Line2D([0], [0], label=r'$\pi^0 (\eta^0) \to \gamma Z^\prime$', color='g')
189194
handles, labels = plt.gca().get_legend_handles_labels()
190195
handles.extend([line_ebrem, line_pbrem, line_pion])
@@ -200,7 +205,7 @@ def plot_muon_limits(mass_ratio=2.1):
200205
plt.xlabel(r"$m_N$ [MeV]", fontsize=14)
201206
plt.xticks(fontsize=14)
202207
plt.yticks(fontsize=14)
203-
plt.xlim((10.0, 1e4))
208+
plt.xlim((50.0, 1e4))
204209
plt.ylim((1e-12, 1.0e-2))
205210
plt.tight_layout()
206211
plt.show()
@@ -228,6 +233,8 @@ def plot_muon_limits_sbn(mass_ratio=2.1):
228233
color='indianred', cls=[3.0*3.0], ls=['dashed'])
229234
c5 = get_Usq_contour_from_file("scans/hnl_DQ_scan_mass-ratio-2_pbrem_umu.txt", color='g', cls=[3.0])
230235
c6 = get_Usq_contour_from_file("scans/hnl_DQ_scan_mass-ratio-2_pbrem_umu.txt", color='g', cls=[3.0e-2], ls=['dashed'])
236+
237+
231238
elif mass_ratio == 5:
232239
c1 = get_Usq_contour_from_file("scans/hnl_SBND_scan_mass-ratio-5_pbrem_umu.txt", color='indianred', cls=[3.0])
233240
c3 = get_Usq_contour_from_file("scans/hnl_MUB_scan_mass-ratio-5_pbrem_umu.txt", color='cadetblue', cls=[3.0])
@@ -236,6 +243,7 @@ def plot_muon_limits_sbn(mass_ratio=2.1):
236243
color='indianred', cls=[3.0*3.0], ls=['dashed'])
237244
c5 = get_Usq_contour_from_file("scans/hnl_DQ_scan_mass-ratio-5_pbrem_umu.txt", color='g', cls=[3.0])
238245
c6 = get_Usq_contour_from_file("scans/hnl_DQ_scan_mass-ratio-5_pbrem_umu.txt", color='g', cls=[3.0e-2], ls=['dashed'])
246+
239247

240248
# plot DQ projections
241249
plt.plot(1e3*dq_18POT_Umu[:,0], dq_18POT_Umu[:,1], linewidth=1.0, color='mediumpurple', ls='dashdot')
@@ -249,10 +257,10 @@ def plot_muon_limits_sbn(mass_ratio=2.1):
249257
plt.text(1750.0, 5e-4, "Laboratory Limits\n(Mixing-only)", fontsize=12)
250258

251259
# text for microboone
252-
plt.text(50, 1.0e-7, "MicroBooNE\n(Mixing Only)", rotation=-45.0, color='teal')
260+
plt.text(50, 1.0e-7, "MicroBooNE\n(Mixing Only)", rotation=-37.0, color='teal')
253261

254262
# text for DQ
255-
plt.text(600.0, 3e-7, "DQ $10^{18}$\n(Mixing Only)", rotation=-25.0, color="mediumpurple")
263+
plt.text(420.0, 5e-7, "DQ $10^{18}$\n(Mixing Only)", rotation=-22.0, color="mediumpurple")
256264
plt.text(2000.0, 7e-9, "DQ $10^{20}$\n(Mixing Only)", rotation=90.0, color="mediumpurple")
257265

258266
line_sbnd = Line2D([0], [0], label=r'SBND', color='indianred')
@@ -263,7 +271,7 @@ def plot_muon_limits_sbn(mass_ratio=2.1):
263271
line_dq_20 = Line2D([0], [0], label=r'DarkQuest ($10^{20}$ POT)', color='g', ls='dashed')
264272
handles, labels = plt.gca().get_legend_handles_labels()
265273
handles.extend([line_sbnd, line_mub, line_icarus, line_sbnd_dump, line_dq, line_dq_20])
266-
plt.legend(handles=handles, loc="lower left", framealpha=1, fontsize=10)
274+
#plt.legend(handles=handles, loc="lower left", framealpha=1, fontsize=10)
267275

268276
plt.yscale('log')
269277
plt.xscale('log')
@@ -275,7 +283,7 @@ def plot_muon_limits_sbn(mass_ratio=2.1):
275283
plt.xlabel(r"$m_N$ [MeV]", fontsize=14)
276284
plt.xticks(fontsize=14)
277285
plt.yticks(fontsize=14)
278-
plt.xlim((10.0, 1e4))
286+
plt.xlim((50.0, 1e4))
279287
plt.ylim((1e-12, 1.0e-2))
280288
plt.tight_layout()
281289
plt.show()
@@ -296,20 +304,30 @@ def plot_tau_limits(mass_ratio=2.1):
296304

297305

298306
if mass_ratio == 2.1:
307+
c4 = get_Usq_contour_from_file("scans/recast_scans/hnl_BEBC_scan_mass-ratio-2_pbrem_utau.txt", color='gray', cls=[3.5])
308+
309+
plt.text(370.0, 8e-5, r"BEBC\\($pp \to X Z^\prime$)", rotation=0.0, color='gray')
310+
299311
c1 = get_Usq_contour_from_file("scans/hnl_DUNE_scan_mass-ratio-2_pbrem_utau.txt", color='indianred')
300312
c2 = get_Usq_upper_lower_limits("scans/hnl_DUNE_scan_mass-ratio-2_ebrem_utau.txt", color='royalblue')
301313
c3 = get_Usq_contour_from_file("scans/hnl_DUNE_scan_mass-ratio-2_meson_utau.txt", color='g')
314+
302315
elif mass_ratio == 5:
303316
c1 = get_Usq_contour_from_file("scans/hnl_DUNE_scan_mass-ratio-5_pbrem_utau.txt", color='indianred')
304317
c2 = get_Usq_upper_lower_limits("scans/hnl_DUNE_scan_mass-ratio-5_ebrem_utau.txt", color='royalblue')
305318
c3 = get_Usq_contour_from_file("scans/hnl_DUNE_scan_mass-ratio-5_meson_utau.txt", color='g')
306319

320+
#c4 = get_Usq_contour_from_file("scans/recast_scans/hnl_BEBC_scan_mass-ratio-5_pbrem_utau.txt", color='k')
321+
322+
323+
# BEBC Label
324+
#plt.text(444.0, 4.4e-7, "BEBC Recast with $Z^\prime$")
307325

308326
# plot seesaw line
309327
masses = np.logspace(0, 4, 100)
310328
plt.fill_between(masses, y1=seesaw_lower(masses), y2=seesaw_upper(masses), color='gold', alpha=0.3)
311329

312-
plt.text(30.0, 1e-9, "Seesaw", rotation=-14.0, fontsize=12)
330+
plt.text(100.0, 4e-10, "Seesaw", rotation=-12.0, fontsize=12)
313331
plt.text(1500.0, 5e-4, "Laboratory Limits\n(Mixing-only)", fontsize=12)
314332

315333
line_pbrem = Line2D([0], [0], label=r'$p N \to p N Z^\prime$', color='indianred')
@@ -330,7 +348,7 @@ def plot_tau_limits(mass_ratio=2.1):
330348
plt.xlabel(r"$m_N$ [MeV]", fontsize=14)
331349
plt.xticks(fontsize=14)
332350
plt.yticks(fontsize=14)
333-
plt.xlim((10.0, 1e4))
351+
plt.xlim((50.0, 1e4))
334352
plt.ylim((1e-12, 1.0e-2))
335353
plt.tight_layout()
336354
plt.show()
@@ -353,6 +371,11 @@ def plot_tau_limits_sbn(mass_ratio=2.1):
353371
color='indianred', cls=[3.0*3.0], ls=['dashed'])
354372
c5 = get_Usq_contour_from_file("scans/hnl_DQ_scan_mass-ratio-2_pbrem_utau.txt", color='g', cls=[3.0])
355373
c6 = get_Usq_contour_from_file("scans/hnl_DQ_scan_mass-ratio-2_pbrem_utau.txt", color='g', cls=[3.0e-2], ls=['dashed'])
374+
375+
c7 = get_Usq_contour_from_file("scans/recast_scans/hnl_BEBC_scan_mass-ratio-2_pbrem_utau.txt", color='gray', cls=[3.5])
376+
377+
plt.text(330.0, 7e-7, r"BEBC\\($pp \to X Z^\prime$)", rotation=0.0, color='gray', fontsize=10)
378+
356379
elif mass_ratio == 5:
357380
c1 = get_Usq_contour_from_file("scans/hnl_SBND_scan_mass-ratio-5_pbrem_utau.txt", color='indianred', cls=[3.0])
358381
c3 = get_Usq_contour_from_file("scans/hnl_MUB_scan_mass-ratio-5_pbrem_utau.txt", color='cadetblue', cls=[3.0])
@@ -362,6 +385,7 @@ def plot_tau_limits_sbn(mass_ratio=2.1):
362385
c5 = get_Usq_contour_from_file("scans/hnl_DQ_scan_mass-ratio-5_pbrem_utau.txt", color='g', cls=[3.0])
363386
c6 = get_Usq_contour_from_file("scans/hnl_DQ_scan_mass-ratio-5_pbrem_utau.txt", color='g', cls=[3.0e-2], ls=['dashed'])
364387

388+
365389
# plot DQ projections
366390
plt.plot(1e3*dq_18POT_Utau[:,0], dq_18POT_Utau[:,1], linewidth=1.0, color='mediumpurple', ls='dashdot')
367391
plt.plot(1e3*dq_20POT_Utau[:,0], dq_20POT_Utau[:,1], linewidth=1.0, color='mediumpurple', ls='dotted')
@@ -375,11 +399,11 @@ def plot_tau_limits_sbn(mass_ratio=2.1):
375399

376400
# text for DQ
377401
if mass_ratio == 2.1:
378-
plt.text(200.0, 1e-4, "DQ $10^{18}$\n(Mixing Only)", rotation=-32.0, color="mediumpurple")
379-
plt.text(200.0, 2.5e-6, "DQ $10^{20}$\n(Mixing Only)", rotation=-32.0, color="mediumpurple")
402+
plt.text(200.0, 1e-4, "DQ $10^{18}$\n(Mixing Only)", rotation=-32.0, color="mediumpurple", zorder=8)
403+
plt.text(200.0, 4e-6, "DQ $10^{20}$\n(Mixing Only)", rotation=-28.0, color="mediumpurple", zorder=8)
380404
else:
381-
plt.text(474.0, 1e-5, "DQ $10^{18}$\n(Mixing Only)", rotation=-30.0, color="mediumpurple")
382-
plt.text(1750.0, 5e-7, "DQ $10^{20}$\n(Mixing Only)", rotation=90.0, color="mediumpurple")
405+
plt.text(474.0, 1e-5, "DQ $10^{18}$\n(Mixing Only)", rotation=-30.0, color="mediumpurple", zorder=8)
406+
plt.text(1750.0, 5e-7, "DQ $10^{20}$\n(Mixing Only)", rotation=90.0, color="mediumpurple", zorder=8)
383407

384408
line_sbnd = Line2D([0], [0], label=r'SBND', color='indianred')
385409
line_sbnd_dump = Line2D([0], [0], label=r'SBND Dump Mode', color='indianred', ls='dashed')
@@ -389,7 +413,9 @@ def plot_tau_limits_sbn(mass_ratio=2.1):
389413
line_dq_20 = Line2D([0], [0], label=r'DarkQuest ($10^{20}$ POT)', color='g', ls='dashed')
390414
handles, labels = plt.gca().get_legend_handles_labels()
391415
handles.extend([line_sbnd, line_mub, line_icarus, line_sbnd_dump, line_dq, line_dq_20])
392-
plt.legend(handles=handles, loc="lower left", framealpha=1, fontsize=10)
416+
417+
if mass_ratio == 5:
418+
plt.legend(handles=handles, loc="lower left", framealpha=1, fontsize=10)
393419

394420
plt.yscale('log')
395421
plt.xscale('log')
@@ -401,17 +427,17 @@ def plot_tau_limits_sbn(mass_ratio=2.1):
401427
plt.xlabel(r"$m_N$ [MeV]", fontsize=14)
402428
plt.xticks(fontsize=14)
403429
plt.yticks(fontsize=14)
404-
plt.xlim((10.0, 1e4))
430+
plt.xlim((50.0, 1e4))
405431
plt.ylim((1e-12, 1.0e-2))
406432
plt.tight_layout()
407433
plt.show()
408434

409435

410436
def main():
411437
plot_muon_limits(mass_ratio=2.1)
412-
plot_tau_limits(mass_ratio=2.1)
413-
plot_muon_limits(mass_ratio=5)
414-
plot_tau_limits(mass_ratio=5)
438+
#plot_tau_limits(mass_ratio=2.1)
439+
#plot_muon_limits(mass_ratio=5)
440+
#plot_tau_limits(mass_ratio=5)
415441

416442
#plot_muon_limits_sbn(mass_ratio=2.1)
417443
#plot_tau_limits_sbn(mass_ratio=2.1)

0 commit comments

Comments
 (0)