In [5]:
import dmdd
import numpy as np
In [ ]:
class Simulation(object):
"""
A simulation of dark-matter direct-detection data under a given experiment and scattering model.
This object handles a single simulated data set (nuclear recoil energy
spectrum). It is generaly initialized and used by the :class:`MultinestRun`
object, but can be used stand-alone.
Simulation data will only be generated if a simulation with the right
parameters and name does not already exist, or if ``force_sim=True`` is
provided upon :class:`Simulation` initialization; if the data exist, it will
just be read in. (Data is a list of nuclear recoil energies of
"observed" events.) Initializing :class:`Simulation` with given parameters
for the first time will produce 3 files, located by default at
``$DMDD_PATH/simulations`` (or ``./simulations`` if ``$DMDD_PATH`` not
defined):
- .dat file with a list of nuclear-recoil energies (keV), drawn from a
Poisson distribution with an expected number of events given by the
underlying scattering model.
- .pkl file with all relevant initialization parameters for record
- .pdf plot of the simulated recoil-energy spectrum with simulated
data points (with Poisson error bars) on top of the underlying model
:param name:
Identifier for simulation (e.g. 'sim1')
:type name: ``str``
:param experiment:
Experiment for simulation.
:type experiment: :class:`Experiment`
:param model:
Model under which to simulate data.
:type model: :class:`Model`
:param parvals:
Values of model parameters. Must contain the
same parameters as ``model``.
:type parvals: ``dict``
:param path:
The path under which to store the simulations.
:type path: ``str``
:param force_sim:
If ``True``, then redo the simulations no matter what. If ``False``,
then the simulations will be redone if and only if the given simulation
parameters don't match what has already been simulated
for this simulation name.
:type force_sim: ``bool``
:param asimov:
Do asimov simulations. Not currently implemented.
:param nbins_asimov:
Number of asimov bins.
:param plot_nbins:
Number of bins to bin data in for rate plot.
:param plot_theory:
Whether to plot the "true" theoretical rate curve along
with the simulated data.
:param silent:
If ``True``, then print messages will be suppressed.
"""
# defines/initializes simulation here
def __init__(self, name, experiment,
model, parvals,
path=SIM_PATH, force_sim=False,
asimov=False, nbins_asimov=20,
plot_nbins=20, plot_theory=True,
silent=False):
#silent suppresses print messages
self.silent = silent
#needs the dictionary of the parameters varying with the model
if not set(parvals.keys())==set(model.param_names):
raise ValueError('Must pass parameter value dictionary corresponding exactly to model.param_names')
self.model = model #underlying model- passed to function
self.experiment = experiment #underlying experiment- passed to function
#build param_values from the dictionary parvals
self.param_values = [parvals[par] for par in model.param_names]
self.param_names = list(self.model.param_names)
for k,v in self.model.fixed_params.items():
self.param_values.append(v)
self.param_names.append(k)
self.name = name
self.path = path
self.asimov = asimov
self.nbins_asimov = nbins_asimov
#how the files are named when they are saved with the data
self.file_basename = '{}_{}'.format(name,experiment.name)
inds = np.argsort(self.param_names)
sorted_parnames = np.array(self.param_names)[inds]
sorted_parvals = np.array(self.param_values)[inds]
for parname, parval in zip(sorted_parnames, sorted_parvals):
self.file_basename += '_{}_{:.2f}'.format(parname, parval)
#calculate total expected rate
dRdQ_params = model.default_rate_parameters.copy()
allpars = model.default_rate_parameters.copy()
allpars['simname'] = self.name
for i,par in enumerate(model.param_names): #model parameters
dRdQ_params[par] = self.param_values[i]
allpars[par] = self.param_values[i]
for kw,val in experiment.parameters.iteritems(): #add experiment parameters
allpars[kw] = val
dRdQ_params['element'] = experiment.element
self.dRdQ_params = dRdQ_params
#need a little more clarification for what is going on here
#just setting up the equations to be calculated later?
self.model_Qgrid = np.linspace(experiment.Qmin,experiment.Qmax,1000)#start, stop, total number
efficiencies = experiment.efficiency(self.model_Qgrid)
self.model_dRdQ = self.model.dRdQ(self.model_Qgrid,**dRdQ_params)
R_integrand = self.model_dRdQ * efficiencies
self.model_R = np.trapz(R_integrand,self.model_Qgrid) # R recoil
self.model_N = self.model_R * experiment.exposure * YEAR_IN_S # N number of EXPECTED events
#create dictionary of all parameters relevant to simulation
self.allpars = allpars
self.allpars['experiment'] = experiment.name
#record dictionary for relevant coupling normalizations only:
norm_dict = {}
for kw in model.param_names:
if kw in PAR_NORMS:
norm_dict[kw] = PAR_NORMS[kw]
self.allpars['norms'] = norm_dict
#data names to be saved
self.datafile = '{}/{}.dat'.format(self.path,self.file_basename)
self.plotfile = '{}/{}.pdf'.format(self.path,self.file_basename)
self.picklefile = '{}/{}.pkl'.format(self.path,self.file_basename)
#control to make sure simulations are forced if they need to be
if os.path.exists(self.picklefile) and os.path.exists(self.datafile):
fin = open(self.picklefile,'rb')
allpars_old = pickle.load(fin)
fin.close()
if not compare_dictionaries(self.allpars,allpars_old):
print('Existing simulation does not match current parameters. Forcing simulation.\n\n')
force_sim = True
else:
print 'Simulation data and/or pickle file does not exist. Forcing simulation.\n\n'
force_sim = True
if force_sim:
if asimov:
raise ValueError('Asimov simulations not yet implemented!')
else:
Q = self.simulate_data()
np.savetxt(self.datafile,Q)
fout = open(self.picklefile,'wb')
pickle.dump(self.allpars,fout)
fout.close()
self.Q = np.atleast_1d(Q)
self.N = len(self.Q)
else:
if asimov:
raise ValueError('Asimov simulations not yet implemented!')
else:
Q = np.loadtxt(self.datafile)
self.Q = np.atleast_1d(Q)
self.N = len(self.Q) #number of events in simulation
if asimov:
raise ValueError('Asimov not yet implemented!')
else:
self.N = len(self.Q) #number of events in simulation
if force_sim or (not os.path.exists(self.plotfile)): #prints the plot and saves it
self.plot_data(plot_nbins=plot_nbins, plot_theory=plot_theory, save_plot=True)
else:
if not self.silent: #prints the number of events and espected events
print "simulation had %i events (expected %.0f)." % (self.N,self.model_N)
def simulate_data(self):
"""
Do Poisson simulation of data according to scattering model's dR/dQ.
"""
Nexpected = self.model_N
if Nexpected > 0:
npts = 10000 #number of points
Nevents = poisson.rvs(Nexpected) #actual number of events generated from poisson dist
Qgrid = np.linspace(self.experiment.Qmin,self.experiment.Qmax,npts)
efficiency = self.experiment.efficiency(Qgrid)
#probability distribution function
pdf = self.model.dRdQ(Qgrid,**self.dRdQ_params) * efficiency / self.model_R
#cumulative distribution function
cdf = pdf.cumsum()
cdf /= cdf.max()
u = random.rand(Nevents) # Nevents number of random numbers
Q = np.zeros(Nevents)
for i in np.arange(Nevents):
Q[i] = Qgrid[np.absolute(cdf - u[i]).argmin()]
else:
Q = np.array([])
Nevents = 0
Nexpected = 0
if not self.silent:
print "simulated: %i events (expected %.0f)." % (Nevents,Nexpected)
return Q
def plot_data(self, plot_nbins=20, plot_theory=True, save_plot=True,
make_plot=True, return_plot_items=False):
"""
Plot simuated data.
:param plot_nbins:
Number of bins for plotting.
:param plot_theory:
Whether to overplot the theory rate curve on top of the
data points.
:param save_plot:
Whether to save plot under ``self.plotfile``.
:param make_plot:
Whether to make the plot. No reason really to ever
be false unless you only want the "plot items"
returned if ``return_plot_items=True`` is passed.
:param return_plot_items:
If ``True``, then function will return lots of things.
"""
Qhist,bins = np.histogram(self.Q,plot_nbins)
Qbins = (bins[1:]+bins[:-1])/2.
binsize = Qbins[1]-Qbins[0] #valid only for uniform gridding.
Qwidths = (bins[1:]-bins[:-1])/2.
xerr = Qwidths
yerr = Qhist**0.5
Qhist_theory = self.model_dRdQ*binsize*self.experiment.exposure*YEAR_IN_S*self.experiment.efficiency(self.model_Qgrid)
Qbins_theory = self.model_Qgrid
if make_plot:
plt.figure()
plt.title('%s (total events = %i)' % (self.experiment.name,self.N), fontsize=18)
xlabel = 'Nuclear recoil energy [keV]'
ylabel = 'Number of events'
ax = plt.gca()
fig = plt.gcf()
xlabel = ax.set_xlabel(xlabel,fontsize=18)
ylabel = ax.set_ylabel(ylabel,fontsize=18)
if plot_theory:
if self.model.name in MODELNAME_TEX.keys():
label='True model ({})'.format(MODELNAME_TEX[self.model.name])
else:
label='True model'
plt.plot(Qbins_theory, Qhist_theory,lw=3,
color='blue',
label=label)
plt.errorbar(Qbins, Qhist,xerr=xerr,yerr=yerr,marker='o',color='black',linestyle='None',label='Simulated data')
plt.legend(prop={'size':20},numpoints=1)
if save_plot:
plt.savefig(self.plotfile, bbox_extra_artists=[xlabel, ylabel], bbox_inches='tight')
if return_plot_items:
return Qbins, Qhist, xerr, yerr, Qbins_theory, Qhist_theory, binsize
In [ ]: