committing changes from long ago - not documented

constants/constants.py

@@ -56,11 +56,35 @@ class tunlSSA_CsI(object):
         csiDiameter = 2.341 # cm, measured by calipers sep8 2014 (picture)
         standoff_TUNLruns = colliToCsI + csiToZero + csiDiameter + tipToColli
 
+    class tunlSSA_CsI_oneBD(object):
+        """Distances associated with the "one-BD CsI QF run at TUNL SSA
+
+        assumes many things are shared from other SSA run(s)
+        """
+        # cm, distance from tipof gas cell to 0deg face
+        cellToZero = 518.055
+        cellLength = 2.86  # cm, length of gas cell
+        zeroDegLength = 3.81  # cm, length of 0deg detector
+        tipToColli = 148.4 # cm, distance from cell tip to collimator exit
+        colliToZero = 233.8 # cm, distance from collimator exit to face of
+            # 0deg AT ITS CLOSEST LOCATION
+        delta1 = 412.3 - 351.3 # cm, difference between close and mid 0degree loc.
+        delta2 = 444.5 - 412.3 # cm, difference between mid and far 0deg loc
+
+        standoffClose = 351.3 
+        standoffMid = standoffClose + delta1
+        standoffFar = standoffMid + delta2
+
+        colliToCsI = 59.45 #from collimator face to closest side of CsI
+        csiToZero = 355.7 # from closest side of CsI to face of zero degree
+        csiDiameter = 2.341 # cm, measured by calipers sep8 2014 (picture)
+        standoff_TUNLruns = colliToCsI + csiToZero + csiDiameter + tipToColli
+
     class cosine_NaI(object):
         """
         Distances associated with COSINE-collab-affiliated NaI QF measurements
         """
-        
+
 
 
 class qValues(object):
@@ -72,7 +96,27 @@ class qValues(object):
 
 class tofWindows(object):
     """Windows of TOF used for different standoff distances
+
+    Really need to clarify documentation here
+
+    Right now, these are specific to 2016 COHERENT CsI[Na] experiments.
+    This should be its own specialized subclass of tofWindows
     """
     nBins = {'close': 45, 'mid': 50, 'far': 70, 'production': 65}
     maxRange = {'close': 175.0, 'mid': 225.0, 'far': 260.0, 'production': 260.0}
-    minRange = {'close':130.0, 'mid': 175.0, 'far': 190.0, 'production':195.0}
+    minRange = {'close':130.0, 'mid': 175.0, 'far': 190.0, 'production':195.0}
+
+    class csi_oneBD(object):
+        """
+        TOF windows specific to CsI[Na] "one BD" experiment
+        """
+        minRange = {'close': 100, 'mid': 120, 'far': 128}
+        maxRange = {'close': 180, 'mid': 200, 'far': 220}
+        # really, bins should be determined based on range.. i think
+        # at least for this experiment, there's pretty coarse binning
+        # should have some error checking here..
+        # but SO RUSHED so just force into ints
+        closeBinning = int((maxRange['close'] - minRange['close']) / 4)
+        midBinning = int((maxRange['mid'] - minRange['mid']) / 4)
+        farBinning = int((maxRange['far'] - minRange['far'] ) / 4)
+        nBins = {'close': closeBinning, 'mid': midBinning, 'far': farBinning}

tests/advIntermediateTOFmodel.py

@@ -53,9 +53,9 @@
 tof_maxRange = tofWindowSettings.maxRange[standoffName[runNumber]]
 tof_range = (tof_minRange,tof_maxRange)
 
-eD_bins = 150
+eD_bins = 240
 eD_minRange = 200.0
-eD_maxRange = 1700.0
+eD_maxRange = 2600.0
 eD_range = (eD_minRange, eD_maxRange)
 eD_binSize = (eD_maxRange - eD_minRange)/eD_bins
 eD_binCenters = np.linspace(eD_minRange + eD_binSize/2,
@@ -65,7 +65,7 @@
 
 x_bins = 100
 x_minRange = 0.0
-x_maxRange = distances.tunlSSA_CsI.cellLength
+x_maxRange = distances.tunlSSA_CsI_oneBD.cellLength
 x_range = (x_minRange,x_maxRange)
 x_binSize = (x_maxRange - x_minRange)/x_bins
 x_binCenters = np.linspace(x_minRange + x_binSize/2,
@@ -78,8 +78,8 @@
 
 
 # PARAMETER BOUNDARIES
-min_e0, max_e0 = 600.0,1100.0
-min_sigma_0,max_sigma_0 = 0.02, 0.3
+min_e0, max_e0 = 1000.0,2600.0
+min_sigma_0,max_sigma_0 = 0.02, 0.5
 
 
 

tests/testPPC.py

@@ -10,6 +10,7 @@
 from utilities.ionStopping import ionStopping
 from utilities.ppcTools import ppcTools
 from scipy.stats import (skewnorm, lognorm)
+from matplotlib import rcParams
 import argparse
 
 
@@ -19,8 +20,8 @@
 
 
 argParser = argparse.ArgumentParser()
-argParser.add_argument('-file')
-argParser.add_argument('-tofDataFile')
+argParser.add_argument('-file', default='/Users/grayson/Documents/mcscratch/burninchain.dat')
+argParser.add_argument('-tofDataFile', default='/Users/grayson/Documents/quenchingFactors/2016-Jan/coda/multistandoff.dat')
 argParser.add_argument('-nParamSamples', default=50, type=int)
 argParser.add_argument('-nBinsE', default=100, type=int)
 argParser.add_argument('-nBinsX', default=20, type=int)
@@ -35,6 +36,10 @@
 
 ppcTool = ppcTools(chainFilename, nSamplesFromTOF=5000,
                    nBins_eD = nBinsE, nBins_x = nBinsX, nRuns = numRuns)
+ppcTool.paramNames = [r'$E_{beam}$ (keV)', r'$\theta$ (keV)', r'$m$ (keV)', 
+                      r'$\sigma$ (1 / keV)', r'$N_1$ (counts)', 
+                      r'$N_2$ (counts)', r'$N_3$ (counts)', r'$N_4$ (counts)',
+                      r'$N_5$ (counts)']
 #ppc, ppcNeutronSpectra = ppcTool.generatePPC( nChainEntries = nParamSamples )
 returnVal = ppcTool.generatePPC( nChainEntries = nParamSamples )
 ppc = returnVal[0]
@@ -102,9 +107,10 @@
 runColors = ["#ff678f", "#004edc", "#e16400", "#1fceff", "#781616", "#149c00"] 
 runNames=['Middle','Close', 'Close (detuned)', 'Far', 'Production']         
 tofXvals = [np.linspace(minT, maxT, bins) for minT, maxT, bins in zip(ppcTool.tof_minRange, ppcTool.tof_maxRange, ppcTool.tofRunBins)]
-fig, axes = plt.subplots(numRuns)
+fig, axes = plt.subplots(numRuns, figsize =(8.5,10.51))
 plotIndexOrder = [1,2,0,3,4]
-ylims=[30e3,40e3,15e3,40e3]
+ylims=[30e3,12e3,22e3,35e3]
+xlims=[(130,175),(130,175),(175,225),(190,260)]
 for idx in range(numRuns):
     index = plotIndexOrder[idx]
     ax=axes[idx]
@@ -122,13 +128,15 @@
      #       statsAllData[plotIndexOrder[idx]][2,:], color='red', alpha=0.4)
     ax.set_ylabel('Counts')
     ax.legend(loc='upper right')
+    ax.set_ylim(0., ylims[idx])
+    ax.set_xlim(xlims[idx])
     #ax.set_ylim(0, ylims[idx])
     #ax.set_xbound(ppcTool.tof_minRange[plotIndexOrder[idx]], 
      #           ppcTool.tof_minRange[plotIndexOrder[idx]])
 #plt.tight_layout()
 plt.xlabel('Time of flight (ns)')
 plt.draw()
-plt.savefig('PPC_on_data.png', dpi=400)
+plt.savefig('PPC_on_data.pdf')
 
 
 
@@ -144,15 +152,15 @@
 dColors = ['#a6611a','#dfc27d']
 nColors = ['#018571','#80cdc1']
 eN_xVals = ppcTool.eN_binCenters
-fig, axNeutrons = plt.subplots()
+fig, axNeutrons = plt.subplots(figsize=(8.5,5.253))
 
 nBand= axNeutrons.fill_between(eN_xVals, neutronStats[0,:], neutronStats[2,:],
                         facecolor=nColors[0], linewidth=0, alpha=0.4,
                         label='Neutron energies')
 axNeutrons.plot(eN_xVals, neutronStats[1,:], color=nColors[0] )
 axNeutrons.tick_params('x', colors=nColors[0])
 axNeutrons.set_xlabel('Neutron energy (keV)', color=nColors[0])
-axNeutrons.set_ylabel('Intensity')
+axNeutrons.set_ylabel('Intensity (a.u.)')
 
 axDeuterons = axNeutrons.twiny()#.twinx()
 
@@ -177,16 +185,18 @@
 #axDeuterons.set_ylabel('Counts (a.u.)')
 features = [dZBand, dBand, nBand]
 labels = [feat.get_label() for feat in features]
-axNeutrons.legend(features, labels, loc=0)
+axNeutrons.legend(features, labels, loc=0,
+                  labelspacing=0.75)
 plt.draw()
-plt.savefig('PPC_neutronSpec.png', dpi=400)
-
+plt.savefig('PPC_neutronSpec.pdf')
 
 
 
 
 
+rcParams.update({'figure.autolayout':False})
 #ppcTool.makeCornerPlot(plotFilename = 'corner_allParams.png')
-#ppcTool.makeCornerPlot(paramIndexHigh = 4, plotFilename = 'corner_eParams.png')
+#cornerFig, cornerAxes = plt.subplots(4,4, figsize=(8,8), )
+ppcTool.makeCornerPlot(paramIndexHigh = 4, nStepsToInclude=150, plotFilename = 'corner_eParams.pdf')
 
 plt.show()

utilities/ionStopping.py

@@ -42,9 +42,16 @@ def __init__(self, params):
                 Z, A, rho, charge, excitation = params
                 # go ahead and calculate number density of the stopping medium
                 # it won't change and will be used for every dEdx calc
-                eNumDensity = (sciConsts.Avogadro * self.protonNumber *
-                              self.density /
-                              (self.massNumber*physics.molarMassConstant))
+
+                # this was originally here
+                # the "self" references don't.. seem right 
+                # left here in case i'm missing something
+                #eNumDensity = (sciConsts.Avogadro * self.protonNumber *
+                #              self.density /
+                #              (self.massNumber*physics.molarMassConstant))
+                eNumDensity = (sciConsts.Avogadro * Z *
+                              rho /
+                              (A * physics.molarMassConstant))
 
                 self.materials.append([Z, A, rho, eNumDensity, excitation])
 

utilities/ppcTools.py

@@ -111,14 +111,16 @@ def __init__(self, chainFilename, nSamplesFromTOF, nBins_eD = 100, nBins_x = 20,
 
 
     def generateModelData(self, params, standoffDistance, range_tof, nBins_tof, ddnXSfxn,
-                      dedxfxn, beamTimer, nSamples, getPDF=False):
+                      dedxfxn, beamTimer, nSamples, getPDF=False, storeDTOF=False):
         """
         Generate model data with cross-section weighting applied
         ddnXSfxn is an instance of the ddnXSinterpolator class -
         dedxfxn is a function used to calculate dEdx -
         probably more efficient to these in rather than reinitializing
         one each time
         This is edited to accommodate multiple standoffs being passed
+        
+        storeDTOF is an option added to add capabilities in MCNP SDEF generation
         """
         beamE, eLoss, scale, s, scaleFactor = params
         e0mean = 900.0
@@ -163,6 +165,7 @@ def generateModelData(self, params, standoffDistance, range_tof, nBins_tof, ddnX
         tofWeights = []
         eN_list = []
         eN_atEachX = np.zeros(self.eD_bins)
+#        dtofs = []
         for index, weight in np.ndenumerate( drawHist2d ):
             cellLocation = self.x_binCenters[index[0]]
             effectiveDenergy = (e0mean + self.eD_binCenters[index[1]])/2
@@ -175,6 +178,8 @@ def generateModelData(self, params, standoffDistance, range_tof, nBins_tof, ddnX
             tofs.append( tof_d + tof_n + zeroD_times )
             tofWeights.append(weight * zeroD_weights)
             eN_list.append(weight)
+#            if storeDTOF:
+#                dtofs.append(tof_d)
             if index[1] == self.eD_binMax:
                 eN_arr = np.array(eN_list)
                 eN_atEachX = np.vstack((eN_atEachX, eN_arr))
@@ -348,6 +353,46 @@ def sampleInitialEnergyDist(self, nSamples = 100, returnNormed=False):
             return dZeroSamples[1:]*self.eD_binSize
         return dZeroSamples[1:]
 
+
+
+    def getDTOFdistribution(self):
+        '''Produce a distribution of deuteron time-of-flight through 
+        gas cell from PPC'''
+        totalChainSamples = len(self.chain[-50:,:,0].flatten())
+
+        dtofHist = np.zeros((self.x_bins, 100))
+
+        dedxForODE = lambda x, y: self.stoppingModel.dEdx(energy=y,x=x)
+
+        samplesToGet = np.random.randint( 0, totalChainSamples, 1 )
+        for sampleToGet in samplesToGet:
+            modelParams = []
+            for nParam in range(self.nParams):
+                modelParams.append(self.chain[-50:,:,nParam].flatten()[sampleToGet])
+
+            e0, loc, scale, s = modelParams[:4]
+
+            eZeros = np.repeat( e0, 1000 )
+            eZeros -= lognorm.rvs(s=s, loc=loc, scale=scale, size=1000)
+            print(eZeros)
+            odesolver = ode( dedxForODE ).set_integrator('dopri5')
+            odesolver.set_initial_value(eZeros)
+            eD_atEachX = np.zeros(self.eD_bins)
+            res = None
+            evPoints = []
+            for idx, xEvalPoint in enumerate(self.x_binCenters):
+                sol = odesolver.integrate(xEvalPoint)
+                print('x loc: {}\t{}'.format(xEvalPoint,sol[0]))
+                if idx == 0:
+                    res = sol
+                else:
+                    res = np.vstack((res, sol))
+                evPoints.append(xEvalPoint)
+        odesolver.set_initial_value([900.0,850.0])
+        testSols = [odesolver.integrate([1.0,1.0]), odesolver.integrate([2.0,2.0])]
+        print('test solutions...\ntest 1: {}\ntest 2: {}'.format(testSols[0], testSols[1]))
+        return res, evPoints, sol
+
 
     def makeSDEF_sia_cumulative(self, distNumber = 100):
         """Produce an MCNP SDEF card for the neutron distribution, marginalized over the cell length
@@ -377,15 +422,15 @@ def makeSDEF_sia_cumulative(self, distNumber = 100):
         return self.sdef_sia_cumulative             
 
 
-    def makeCornerPlot(self, paramIndexLow = 0, paramIndexHigh = None, plotFilename = 'ppcCornerOut.png'):
+    def makeCornerPlot(self, paramIndexLow = 0, paramIndexHigh = None, nStepsToInclude = 50, plotFilename = 'ppcCornerOut.png'):
         """Produce a corner plot of the parameters from the chain
         
         paramIndexLow and paramIndexHigh define the upper and lower boundaries of the parameter indices to plot"""
         if paramIndexHigh == None:
             paramIndexHigh = self.nParams
-        samples = self.chain[-50:,:,paramIndexLow:paramIndexHigh].reshape((-1, paramIndexHigh - paramIndexLow))
+        samples = self.chain[-nStepsToInclude:,:,paramIndexLow:paramIndexHigh].reshape((-1, paramIndexHigh - paramIndexLow))
         import corner as corn
         cornerFig = corn.corner(samples, labels=self.paramNames[paramIndexLow:paramIndexHigh], 
-                                quantiles=[0.15, 0.5, 0.84], show_titles=True,
+                                quantiles=[0.16, 0.5, 0.84], show_titles=True,
                                 title_kwargs={'fontsize': 12})
-        cornerFig.savefig(plotFilename, dpi=300)
+        cornerFig.savefig(plotFilename)