Initial commit of code examples for data release

Code accompanying the 2018 release of data from the COHERENT Collaboration

This release corresponds to the result published in arXiv:1708.01294[nucl-ex]

For a full description, see the PDF document included within the release. The release can be found at http://coherent.ornl.gov/data or on zenodo.org, DOI: 10.5281/zenodo.1228631

This code is intended to provide examples of reading the data comprising the release

Included Python scripts are:

	coherent_readDataExample.py: this reads one of the 2-D data files and produces an image of a corresponding 2-D histogram. A function included in this script can be reused to read all 2-D data included.

	coherent_readPromptPDF.py: this reads the text file describing the distribution in photoelectron space of events from prompt, SNS-produced neutrons, producing an image of this distribution with and without analysis efficiency applied.

	coherent_readTiming.py: this reads the text file describing the time distributions for prompt neutrons and both the prompt and delayed CEvNS populations. An image is produced which shows these distributions.

	coherent_readParameters.py: this reads the YAML file which includes experimental parameters or details that can be described by single values with (or without) uncertainties.

#

# coherent_readDataExample.py

#

# This file demonstrates reading of the 2-dimensional data packaged in the 

# COHERENT Collaboration data release associated with arXiv:1708.01294 [nucl-ex]

#

# It will read the beam-on, coincidence-region data and produce a PDF file

# of 2-D histogram of this data, showing also the projections onto both 

# the photoelectron and arrival time axes

# 

#

# Command line arguments specify the data file to read and the image file to produce

# default values allow one to simply run 'python coherent_readDataExample.py' from within the 'code' directory of the release

# this will read data stored in the 'data' directory of the release and produce a test image in 'code'

#

# created 2018 March, Grayson C. Rich

# gcrich@uchicago.edu

#

# 

from __future__ import print_function

import numpy as np

import matplotlib.pyplot as plt

from matplotlib import gridspec

import matplotlib.ticker as ticker

import argparse

def parseCOHERENTdataFile(filename):

    '''

    Read in a data file from the first COHERENT data release, specified by argument filename

    Returns a 2-D numpy array representing the histogrammed data and two 1-D numpy arrays with the bin centers for each axis

    '''

    csvFileData = np.genfromtxt(filename,delimiter=', ')

    # remove duplicate entries and sort

    BinCentersX = sorted(list(set(csvFileData[:,0])))

    BinCentersY = sorted(list(set(csvFileData[:,1])))

    nBinsX = len(BinCentersX)

    nBinsY = len(BinCentersY)

    newHist = np.zeros((nBinsX,nBinsY))

    for yIndex in range(nBinsY):

        for xIndex in range(nBinsX):

            newHist[xIndex,yIndex] = csvFileData[xIndex * nBinsY + yIndex,2]

    return newHist, BinCentersX, BinCentersY

def roundToMultiple(x, base = 5):

    '''

    Round a number to nearest multiple of base

    From: https://stackoverflow.com/questions/2272149/round-to-5-or-other-number-in-python

    '''

    return int(base * round(float(x)/base))

if __name__=='__main__':

        #Parse command line arguments

        argParser = argparse.ArgumentParser()

        argParser.add_argument('-inputFilename', type=str, default='../data/data_coincidence_beamOn.txt',

                                help='Name (with path, if in other directory) of the text file containing the 2-D beam-on, coincidence-region data')

        argParser.add_argument('-outputFilename', type=str, default='./coincidence_beamOn_2D.pdf',

                                help='Name of the file to be produced showing the 2-D histogram of the beam-on, coincidence-region data')

        parsedArgs = argParser.parse_args()

        inputFilename = parsedArgs.inputFilename

        outputFilename = parsedArgs.outputFilename

        #Load the data

        theHist, binCentersX, binCentersY = parseCOHERENTdataFile(inputFilename)

        peBinCenters = binCentersX

        peBinWidth = peBinCenters[1] - peBinCenters[0]

        timeBinCenters = binCentersY

        timeBinWidth = timeBinCenters[1] - timeBinCenters[0]

        #Format the display

        axisLabelSets= [['Photoelectrons','Counts ( / 2 PE )'],

                        [r'Arrival time ($\mu$s)',r'Counts ( / 0.5 $\mu$s )']]

        layout = gridspec.GridSpec(2,2, width_ratios=[3,1], height_ratios=[3,1],

                                   left=0.15, right=0.95, top=0.85, bottom=0.1,

                                   hspace=0.07, wspace=0.07)

        fig = plt.figure(figsize=(6,6))

        ax = fig.add_subplot(layout[0])

        hist2dPlot = ax.imshow(theHist, interpolation='nearest',

                                 cmap='cubehelix_r', aspect='auto',

                                 extent=(0., 12, 50, 0.), origin='upper')

        cbaxis = fig.add_axes([0.15, 0.87, 0.58, 0.03], frameon=False)

        cbar = fig.colorbar(hist2dPlot, cax=cbaxis, orientation='horizontal')

        cbar.ax.xaxis.set_ticks_position('top')

        cbar.ax.xaxis.set_label_position('top')

        cbar.ax.set_xlabel(r'Counts ( / PE$\cdot\mu$s )')

        ax.set_xlabel(r'Arrival time ( / 0.5 $\mu$s)')

        ax.set_ylabel('Photoelectrons (PEs)' )

        ax_pe = fig.add_subplot(layout[1], sharey=ax)

        ax_time = fig.add_subplot(layout[2], sharex=ax)

        ax_pe.bar(left=None, height=peBinWidth, bottom=peBinCenters-peBinWidth/2,

                  width=np.sum(theHist, axis=1), orientation=u'horizontal',

                  color='black', linewidth=0, alpha=0.3)

        ax_time.bar(timeBinCenters-timeBinWidth/2, np.sum(theHist, axis=0),

                    width=timeBinWidth, color='black', linewidth=0, alpha=0.3)

        ax_pe.set_ylim(50.,0.)

        ax_time.set_xlim(0., 12)

        #

        # determine axis ticker label spacing for the projections

        # note that this isn't elegantly refined

        # it isn't intended to work anywhere but here

        #

        # determine ticker widths for PE axis

        maxValue_PEaxis = np.max(np.sum(theHist, axis=1))

        peAxisTickerWidth = roundToMultiple(maxValue_PEaxis/5 + 10, 10)

        peAxisMinorTickerWidth = round(peAxisTickerWidth/5)

        # determine ticker widths for time axis

        maxValue_timeAxis = np.max(np.sum(theHist, axis=0))

        timeAxisTickerWidth = roundToMultiple(maxValue_timeAxis/5 + 10, 10) # alternatively, set to a number of your choosing

        timeAxisMinorTickerWidth = round(timeAxisTickerWidth/5) # same

        ax_time.set_xlabel(axisLabelSets[1][0])

        ax_time.set_ylabel(axisLabelSets[1][1])

        ax_time.yaxis.set_major_locator(ticker.MultipleLocator(timeAxisTickerWidth))

        ax_time.yaxis.set_minor_locator(ticker.MultipleLocator(timeAxisMinorTickerWidth))

        ax_pe.set_xlabel(axisLabelSets[0][1])

        ax_pe.xaxis.set_label_position('top')

        ax_pe.xaxis.tick_top()

        ax_pe.xaxis.set_ticks_position('both')

        ax_pe.xaxis.set_major_locator(ticker.MultipleLocator(peAxisTickerWidth))

        ax_pe.xaxis.set_minor_locator(ticker.MultipleLocator(peAxisMinorTickerWidth))

        ax_pe.xaxis.label.set_size(12)

        ax_time.yaxis.label.set_size(12)

        ax.xaxis.set_minor_locator(ticker.MultipleLocator(1))

        ax.yaxis.set_minor_locator(ticker.MultipleLocator(2))

        plt.setp(ax.get_xticklabels(), visible=False)

        plt.setp(ax_pe.get_yticklabels(), visible=False)

        ax.xaxis.label.set_visible(False)

        plt.savefig(outputFilename)

        plt.show()

#

# coherent_readParameters.py

#

#

# A simple demonstration of reading the YAML parameter file associated 

# with the data release from the COHERENT Collaboration pertaining to arXiv:1708.01294 [nucl-ex]

#

# Command line arguments specify the data file to read

# default values allow one to simply run 'python coherent_readDataExample.py' from within the 'code' directory of the release

# this will read data stored in the 'data' directory of the release

# 

# created 2018 March, Grayson C. Rich

# gcrich@uchicago.edu

#

#

from __future__ import print_function

import yaml

import argparse

argParser = argparse.ArgumentParser()

argParser.add_argument('-inputFilename', type=str, default='../data/coherent_parameters.yaml',

                        help='Name (with path, if in other directory) of the YAML file containing COHERENT experiment parameters')

parsedArgs = argParser.parse_args()

inputFilename = parsedArgs.inputFilename

parameterDictionary = yaml.load(open(inputFilename,'r'))

print('Printing names of the parameters included in YAML file..\n')

for entry in parameterDictionary:

    paramEntry = parameterDictionary.get(entry)

    paramName = paramEntry.get('name')

    paramValue = paramEntry.get('value')

    subparameters = paramEntry.get('parameters')

    if paramValue != None:

        print('{} has value {}'.format(paramName, paramValue))

    if subparameters:

        print('{} has list of parameters!'.format(paramName))

        for subparam in subparameters:

            subparamName = subparam.get('name')

            subparamValue = subparam.get('value')

            print('\t{} has value {}'.format(subparamName, subparamValue))

    print('\n')

 No newline at end of file

#

# coherent_readPromptPDF.py

#

# Intended to provide simple example of reading the prompt neutron PDF

# which is included in the COHERENT Collaboration data release associated 

# with arXiv:1708.01294 [nucl-ex]

#

# Command line arguments specify the data file to read and the image file to produce

# default values allow one to simply run 'python coherent_readPromptPDF.py' from within the 'code' directory of the release

# this will read data stored in the 'data' directory of the release and produce a test image in 'code'

#

# created 2018 March, Grayson C. Rich

# gcrich@uchicago.edu

#

#

import numpy as np

import matplotlib.pyplot as plt

import argparse

argParser = argparse.ArgumentParser()

argParser.add_argument('-inputFilename', type=str, default='../data/promptPDF.txt',

                        help='Name (with path, if in other directory) of the text file containing the prompt neutron PDF')

argParser.add_argument('-outputFilename', type=str, default='./promptNeutronPDF.pdf',

                        help='Name of the file to be produced showing the expected contribution in PE space of prompt neutrons')

parsedArgs = argParser.parse_args()

inputFilename = parsedArgs.inputFilename

outputFilename = parsedArgs.outputFilename

def parse1DpdfFile(filename):

    '''

    Simple function to read in the PDF data for a 1D PDF as formatted for the COHERENT data release

    Reads data from file (name specified as argument) and returns two numpy arrays containing the bin centers and the bin contents

    '''

    csvData = np.genfromtxt(filename,delimiter=',')

    return csvData[:,0],csvData[:,1]

promptBinCenters, promptContents = parse1DpdfFile(inputFilename)

# now we apply the efficiency curve

#  but first we have to work it out!

def csiEfficiency(x, a=0.6655, k=0.4942, shift=10.8507):

    '''

        Evaluate the acceptance efficiency for the CsI[Na] analysis for a signal with a specified number of photoelectrons

        Note that this should be applied to data/models binned with 1PE granularity - the efficiency for a bin is based on this function's value at the bin center

        So, this would be evaluted at, e.g., 0.5 for efficiency of 1PE signals, 1.5 for efficiency of 2PE signals, etc

        '''

    if x < 5:

        return 0.

    acceptance = a / (1.0 + np.exp(-k*(x-shift)))

    if x >= 5 and x < 6:

        acceptance = acceptance * 0.5

    return acceptance

# get an array of the acceptance efficiencies for each bin

binEfficiencies = np.array([csiEfficiency(bincenter) for bincenter in np.linspace(0.5, 49.5, 50)])

# and the make an array of 'accepted' prompt counts, having applied efficiency

efficiencyAppliedPrompts = promptContents * binEfficiencies

# now rebin the distributions to be consistent with the 2-PE bin widths of

# the data included in this release

promptContents_2PE = []

efficiencyAppliedPrompts_2PE = []

for index in range(int(len(promptContents)/2)):

    promptContents_2PE.append( promptContents[index*2] +

                              promptContents[index*2+1] )

    efficiencyAppliedPrompts_2PE.append( efficiencyAppliedPrompts[2*index] +

                                       efficiencyAppliedPrompts[2*index+1] )

# set up the bin centers and bin widths

# we know we have 2-PE bins, and bin centers at 1 PE, 3 PE, 5 PE, etc

promptBinCenters = np.linspace( 1, 49, 25 )

binWidth = 2.

# plot raw prompt distribution with the efficiency-cut distribution overlaid

fig,ax = plt.subplots(figsize=(5.,3.09))

ax.bar( promptBinCenters-binWidth/2, height = promptContents_2PE, 

       width = binWidth, linewidth = 0, label='Prompt neutrons')

ax.bar( promptBinCenters-binWidth/2, height = efficiencyAppliedPrompts_2PE,

       width=binWidth, linewidth=0, color='gray', 

       label='Analysis efficiency applied')

ax.set_xlabel('Photoelectrons')

ax.set_ylabel('Counts / bin / GW$\cdot$hr')

ax.legend(loc='upper right', frameon=False, numpoints=1, labelspacing=0.05,

          handleheight=2.0)

plt.savefig(outputFilename)

plt.draw()

plt.show()

#

# coherent_readTiming.py

#

# Intended to provide simple example of reading the timing distributions for the signal components

# which is included in the COHERENT Collaboration data release associated 

# with arXiv:1708.01294 [nucl-ex]

#

# Command line arguments specify the data file to read and the image file to produce

# default values allow one to simply run 'python coherent_readTiming.py' from within the 'code' directory of the release

# this will read data stored in the 'data' directory of the release and produce a test image in 'code'

#

# created 2018 March, Grayson C. Rich

# gcrich@uchicago.edu

#

#

import numpy as np

import csv

import matplotlib.pyplot as plt

import argparse

argParser = argparse.ArgumentParser()

argParser.add_argument('-inputDirectory', type=str, default='../data/',

                        help='Path to the directory holding the txt files with timing distribution information')

argParser.add_argument('-outputFilename', type=str, default='./exampleTiming.pdf',

                        help='Name of the file to be produced showing the normalized timing distributions for prompt neutrons and CEvNS')

parsedArgs = argParser.parse_args()

inputDirectory = parsedArgs.inputDirectory

outputFilename = parsedArgs.outputFilename

def parse1DpdfFile(filename):

    '''

    Simple function to read in the PDF data for a 1D PDF as formatted for the COHERENT data release

    Reads data from file (name specified as argument) and returns two numpy arrays containing the bin centers and the bin contents

    '''

    binCenters = []

    binContents = []

    with open(filename,'r') as csvFile:

        csvReader = csv.reader(csvFile, delimiter=',', skipinitialspace=True)

        for row in csvReader:

            if row[0].startswith('#'):

                continue

            binCenters.append(float(row[0]))

            binContents.append(float(row[1]))

    return np.array(binCenters), np.array(binContents)

filenameList = ['arrivalTimePDF_promptNeutrons.txt',

                'arrivalTimePDF_promptNeutrinos.txt', 

                'arrivalTimePDF_delayedNeutrinos.txt']

inputFilenames = [inputDirectory + filename for filename in filenameList]

# read in the PDFs for each of the components

# this isn't very pythonic, but it is clear

# note that all bin centers are the same, so the reuse of timingBinCenters

# is ok

timingBinCenters, promptNeutronPDF = parse1DpdfFile(inputFilenames[0])

timingBinCenters, promptNeutrinoPDF = parse1DpdfFile(inputFilenames[1])

timingBinCenters, delayedNeutrinoPDF = parse1DpdfFile(inputFilenames[2])

# need to figure out bin width for plotting

# assume uniform bin width (which we know to be the case)

binWidth = timingBinCenters[1] - timingBinCenters[0]

# we also need the low edge of the bin, so calculate those

binLowEdges = [binCenter - 0.5 * binWidth for binCenter in timingBinCenters]

# we'll plot the neutron timing separate from the neutrino timing

# prompt neutrons will go on top, neutrinos on the bottom

fig, [topAx, bottomAx] = plt.subplots(2, sharex=True)

topAx.bar(binLowEdges, promptNeutronPDF, binWidth, linewidth=0, 

          color='#b66500', label='Prompt neutrons')

topAx.set_ylim(0., 1.)

topAx.set_ylabel('Intensity / bin')

topAx.legend(loc='upper right', frameon=False, numpoints=1,

             labelspacing=0.05, handleheight=2.0)

bottomAx.bar(binLowEdges, promptNeutrinoPDF, binWidth, linewidth=0, 

             color= '#03e19e', label=r'Prompt CE$\nu$NS')

bottomAx.bar(binLowEdges, delayedNeutrinoPDF, binWidth, linewidth=0,

             color='#ff3072', label=r'Delayed CE$\nu$NS', alpha=0.5)

bottomAx.set_ylim(0.,1.)

bottomAx.set_ylabel('Intensity / bin')

bottomAx.legend(loc='upper right', frameon=False, numpoints=1, 

                labelspacing=0.05, handleheight=2.0)

bottomAx.set_xlim(0., 6)

bottomAx.set_xlabel(r'Arrival time ($\mu$s)')

plt.savefig(outputFilename)

plt.draw()

plt.show()

 No newline at end of file