diff options
| -rw-r--r-- | utils/Makefile | 1 | ||||
| -rwxr-xr-x | utils/gen-dark-matter | 342 |
2 files changed, 343 insertions, 0 deletions
diff --git a/utils/Makefile b/utils/Makefile index 348c4f2..ef21a09 100644 --- a/utils/Makefile +++ b/utils/Makefile @@ -12,5 +12,6 @@ install: $(INSTALL) plot-fit-results $(INSTALL_BIN) $(INSTALL) plot-likelihood $(INSTALL_BIN) $(INSTALL) submit-grid-jobs $(INSTALL_BIN) + $(INSTALL) gen-dark-matter $(INSTALL_BIN) .PHONY: install diff --git a/utils/gen-dark-matter b/utils/gen-dark-matter new file mode 100755 index 0000000..007dde3 --- /dev/null +++ b/utils/gen-dark-matter @@ -0,0 +1,342 @@ +#!/usr/bin/env python +# Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago> +# +# This program is free software: you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) +# any later version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for +# more details. +# +# You should have received a copy of the GNU General Public License along with +# this program. If not, see <https://www.gnu.org/licenses/>. +""" +Script to produce MCPL files for simulating self destructing dark matter. The +arguments to the script are the dark matter mediator mass and energy and the +particle IDs of the two decay products. For example, to simulate a dark matter +mediator with a mass of 100 MeV and total energy 1 GeV decaying to an electron +and a positron: + + $ ./gen-dark-matter -M 100.0 -E 1000.0 -p1 20 -p2 21 -o output + +which will create the directory "output_[UUID]" and produce three files: + + 1. A MCPL file output.mcpl + 2. A MCPL header file called mcpl_header.dat + 3. A SNOMAN command file which can be used to simulate the dark matter + +To simulate it with SNOMAN you can run: + + $ cd output_[UUID] + $ snoman.exe -c output.cmd +""" + +from __future__ import print_function, division +import numpy as np +import string +import uuid + +SNOMAN_MASS = { + 20: 0.511, + 21: 0.511, + 22: 105.658, + 23: 105.658 +} + +PSUP_RADIUS = 840.0 # cm + +# generate a UUID to append to all the filenames so that if we run the same job +# twice we don't overwrite the first job +ID = uuid.uuid1() + +SNOMAN_TEMPLATE = \ +""" +$processor_list 'MCO UCL PRU PCK OUT END' + +* set up output file +$output_format $full_ds +$zdab_option $zdab_max_mc +file out 1 @output +* Don't know if this next line is necessary +file MCO 1 ./MC_Atm_Nu_No_Osc_Snoman_Genie10000_X_0_rseed.dat checkpoint=10 +file mco 2 @mcpl +$mc_event_rate -0.003189 $per_sec + + +$prune_mc $keep +$prune_mcpm $keep +$prune_mcvx_source $keep +$prune_mcvx_boundary $drop +$prune_mcvx_interaction $drop +$prune_mcvx_sink $drop +$prune_mcvx_pre_source $drop + +$mcrun 10000 + +* simulate run conditions for run 10000 +$mc_gen_run_cond $on + +* Don't think this is necessary since it's reset in run_mc_atmospherics +$num_events @num_events + +* titles files for run 10000 +titles /sno/mcprod/dqxx/DQXX_0000010000.dat +titles /sno/mcprod/anxx/titles/pca/anxx_pca_0000010623_p2.dat +titles /sno/output/anxx/titles/neutrino/10000-10999/anxx_nu_0000010000_p12.dat + +* set the average number of noise hits per event +* this comes from the autosno generated MC_Atm_Nu_No_Osc_Snoman_Genie10000_X_1_run_mc_atm_nu_genie.cmd file +set bank TRSP 1 word 2 to 2.051309 + +* From activate_atmospherics.cmd + +$killvx 7 +$killvx_neutron 5 +$egs4_ds $off +$store_full_limit 10 +$max_cer_ge_errors 2000 + +* Enable hadron propagation + +titles sno_hadron_list.dat +titles chetc_sno.dat +titles flukaaf_sno.dat +$enable_hadrons $on + +* Load information for muons + +titles music_sno_info.dat +titles music_double_diff_rock.dat +titles muon.dat +titles muon_param.dat +titles photo_dis.dat +$enable_music_calc $off + +@load_d2o_settings.cmd +@run_mc_atmospherics +""".strip() + +MCPL_HEADER=""" +*DO MCPL 1 -i(30I 2I / 2I 17F) -n2000000 +#. +#. This bank is used to run particles through SNOMAN (with file MCO 2 ... command) +#. This file was derived from a .dat file of Christian Nally. +#. +#. Contact: D. Waller (Carleton) +#. +#. Standard Database Header +#. +19750101 0 20380517 03331900 #. 1..4 Intrinsic validity + 0 0 0 #. 5..7 Data type, Task type, Format no. + 0 0 0 #. 8..10 Creation Date, Time, Source Id. +19750101 0 20380517 03331900 #. 11..14 Effective validity + 0 0 #. 15..16 Entry Date Time +4*0 #. 17..20 Spare +1000000*0 #. 21..30 Temporary data (not in database) +#. +#. End of Standard Database Header +#. +#. User Data. +""" + +class MyTemplate(string.Template): + delimiter = '@' + +def rand_sphere2(R): + """ + Generates a random point inside a sphere of radius R. + """ + while True: + pos = np.random.rand(3)*R + + if np.linalg.norm(pos) < R: + break + + return pos + +def rand_sphere(): + """ + Generates a random point on the unit sphere. + """ + u = np.random.rand() + v = np.random.rand() + + phi = 2*np.pi*u + theta = np.arccos(2*v-1) + + dir = np.empty(3,dtype=float) + + dir[0] = np.sin(theta)*np.cos(phi) + dir[1] = np.sin(theta)*np.sin(phi) + dir[2] = np.cos(theta) + + return dir + +def lorentz_boost(x,n,beta): + """ + Performs a Lorentz boost on the 4-vector `x` along the direction `n` at a + velocity `beta` (in units of the speed of light). + + Formula comes from 46.1 in the PDG "Kinematics" Review. See + http://pdg.lbl.gov/2014/reviews/rpp2014-rev-kinematics.pdf. + """ + n = np.asarray(n,dtype=float) + x = np.asarray(x,dtype=float) + + gamma = 1/np.sqrt(1-beta**2) + + # normalize direction + n /= np.linalg.norm(n) + + # get magnitude of `x` parallel with the boost direction + x_parallel = np.dot(x[1:],n) + + # compute the components of `x` perpendicular + x_perpendicular = x[1:] - n*x_parallel + + E = x[0] + E_new = gamma*E - gamma*beta*x_parallel + x_new = -gamma*beta*E + gamma*x_parallel + + x_new = [E_new] + (x_new*n + x_perpendicular).tolist() + + return np.array(x_new) + +def gen_decay(M, E, m1, m2): + """ + Generator yielding a tuple (v1,v2) of 4-vectors for the daughter particles + of a 2-body decay from a massive particle M with total energy E in the lab + frame. + + Arguments: + + M - mass of parent particle (MeV) + E - Total energy of the parent particle (MeV) + m1 - mass of first decay product + m2 - mass of second decay product + """ + while True: + # direction of 1st particle in mediator decay frame + v = rand_sphere() + # calculate momentum of each daughter particle + # FIXME: need to double check my math + p1 = np.sqrt(((M**2-m1**2-m2**2)**2-4*m1**2*m2**2)/(4*M**2)) + E1 = np.sqrt(m1**2 + p1**2) + E2 = np.sqrt(m2**2 + p1**2) + v1 = [E1] + (p1*v).tolist() + v2 = [E2] + (-p1*v).tolist() + # random direction of mediator in lab frame + n = rand_sphere() + beta = np.sqrt(E**2-M**2)/E + v1_new = lorentz_boost(v1,n,beta) + v2_new = lorentz_boost(v2,n,beta) + + yield v1_new, v2_new + +def test_lorentz_boost(): + """ + Super simple test of the lorentz_boost() function to make sure that if we + boost forward by a speed equal to the particle's original speed, it's + 4-vector looks like [mass,0,0,0]. + """ + m1 = 100.0 + beta = 0.5 + p1 = m1*beta/np.sqrt(1-beta**2) + p1 = [np.sqrt(m1**2 + p1**2)] + [p1,0,0] + + p1_new = lorentz_boost(p1,[1,0,0],0.5) + + assert np.allclose(p1_new,[m1,0,0,0]) + +if __name__ == '__main__': + import argparse + from itertools import islice + import sys + import os + + parser = argparse.ArgumentParser("generate MCPL files for self destructing dark matter") + parser.add_argument("-M", type=float, default=100.0, + help="mass of mediator") + parser.add_argument("-E", type=float, default=100.0, + help="total energy of mediator") + parser.add_argument("-p1", type=int, default=20, + help="SNOMAN particle ID for 1st decay product") + parser.add_argument("-p2", type=int, default=21, + help="SNOMAN particle ID for 2nd decay product") + parser.add_argument("-n", type=int, default=100, + help="number of events to generate") + parser.add_argument("-o", "--output", type=str, required=True, + help="output prefix") + args = parser.parse_args() + + if args.p1 not in SNOMAN_MASS: + print("%i is not a valid particle ID" % args.p1,file=sys.stderr) + sys.exit(1) + + m1 = SNOMAN_MASS[args.p1] + + if args.p2 not in SNOMAN_MASS: + print("%i is not a valid particle ID" % args.p2,file=sys.stderr) + sys.exit(1) + + m2 = SNOMAN_MASS[args.p2] + + if args.M < m1 + m2: + print("mediator mass must be greater than sum of decay product masses",file=sys.stderr) + sys.exit(1) + + if args.E < args.M: + print("mediator energy must be greater than or equal to the mass",file=sys.stderr) + sys.exit(1) + + if args.n < 0: + print("number of events must be positive",file=sys.stderr) + sys.exit(1) + + # The format for the MCPL files as read in by SNOMAN is: + # + # Word Type Description + # ---- ---- -------------------------------- + # 1 I Number of events in the list. + # 2 I Number of words per track. (This may vary because + # polarisation is included only for Cerenkov photons). + # j+1 I Number of particles in this event. + # j+2 I Particle type of this particle. + # j+3 F,F,F X,Y,Z of particle. + # j+6 F Time of particle. + # j+7 F Energy of particle. + # j+8 F,F,F U,V,W of particle. + # j+11 F,F,F x,y,z components of the particle's polarisation. + + + new_dir = "%s_%s" % (args.output,ID.hex) + + os.mkdir(new_dir) + os.chdir(new_dir) + + mcpl_filename = args.output + ".mcpl" + + with open(mcpl_filename, "w") as f: + f.write("%i %i\n" % (args.n, 10)) + + for v1, v2 in islice(gen_decay(args.M,args.E,m1,m2),args.n): + pos = rand_sphere2(PSUP_RADIUS) + p1 = np.linalg.norm(v1[1:]) + p2 = np.linalg.norm(v2[1:]) + f.write(" %i %i %f %f %f %f %f %f %f %f\n" % (2,args.p1,pos[0], pos[1], pos[2],0.0, v1[0], v1[1]/p1, v1[2]/p1, v1[3]/p1)) + f.write(" %i %i %f %f %f %f %f %f %f %f\n" % (2,args.p2,pos[0], pos[1], pos[2],0.0, v2[0], v2[1]/p2, v2[2]/p2, v2[3]/p2)) + + # Write out mcpl_header.dat which is needed by the cmd file + with open("mcpl_header.dat", "w") as f: + f.write(MCPL_HEADER) + + template = MyTemplate(SNOMAN_TEMPLATE) + + mcds_filename = args.output + ".mcds" + cmd_filename = args.output + ".cmd" + + with open(cmd_filename, "w") as f: + f.write(template.safe_substitute(output=mcds_filename, mcpl=mcpl_filename, num_events=str(args.n))) |