diff options
Diffstat (limited to 'utils/dm-search')
| -rwxr-xr-x | utils/dm-search | 28 |
1 files changed, 21 insertions, 7 deletions
diff --git a/utils/dm-search b/utils/dm-search index a486593..25fa185 100755 --- a/utils/dm-search +++ b/utils/dm-search @@ -35,7 +35,7 @@ from sddm import printoptions from sddm.utils import fast_cdf, correct_energy_bias from scipy.integrate import quad from sddm.dm import * -from sddm import SNOMAN_MASS +from sddm import SNOMAN_MASS, AV_RADIUS import nlopt # Likelihood Fit Parameters @@ -494,6 +494,18 @@ if __name__ == '__main__': evs.append(get_events(df.filename.values, merge_fits=True, nhit_thresh=args.nhit_thresh)) ev = pd.concat(evs) + livetime = 0.0 + livetime_pulse_gt = 0.0 + for ev in evs: + if not np.isnan(ev.attrs['time_10_mhz']): + livetime += ev.attrs['time_10_mhz'] + else: + livetime += ev.attrs['time_pulse_gt'] + livetime_pulse_gt += ev.attrs['time_pulse_gt'] + + print("livetime = %.2f" % livetime) + print("livetime (pulse gt) = %.2f" % livetime_pulse_gt) + ev = correct_energy_bias(ev) # Note: We loop over the MC filenames here instead of just passing the @@ -530,10 +542,12 @@ if __name__ == '__main__': ev_mc = ev_mc[~np.isnan(ev_mc.fmin)] muon_mc = muon_mc[~np.isnan(muon_mc.fmin)] - # require (r/r_psup)^3 < 0.9 - ev = ev[ev.r_psup < 0.9] - ev_mc = ev_mc[ev_mc.r_psup < 0.9] - muon_mc = muon_mc[muon_mc.r_psup < 0.9] + # require (r < av radius) + ev = ev[ev.r < AV_RADIUS] + ev_mc = ev_mc[ev_mc.r < AV_RADIUS] + muon_mc = muon_mc[muon_mc.r < AV_RADIUS] + + fiducial_volume = (4/3)*np.pi*(AV_RADIUS)**3 # require psi < 6 ev = ev[ev.psi < 6] @@ -702,11 +716,11 @@ if __name__ == '__main__': fig = plt.figure() for color, dm_particle_id in zip(('C0','C1'),(2020,2222)): - plt.plot(DM_MASSES[dm_particle_id],limits[dm_particle_id],color=color,label='$' + ''.join([particle_id[int(''.join(x))] for x in grouper(str(dm_particle_id),2)]) + '$') + plt.plot(DM_MASSES[dm_particle_id],limits[dm_particle_id]/fiducial_volume/livetime,color=color,label='$' + ''.join([particle_id[int(''.join(x))] for x in grouper(str(dm_particle_id),2)]) + '$') plt.gca().set_xscale("log") despine(fig,trim=True) plt.xlabel("Energy (MeV)") - plt.ylabel("Event Rate Limit (events)") + plt.ylabel("Event Rate Limit (events/cm^3/s)") plt.legend() plt.tight_layout() |