Feature/dark news table (#37)

* reorganizing examples and earthparams files

* starting to add fill interpolation method

* Flush out fill table methods for LIDarkNews

* Minor cleanup

python/LIDarkNews.py

@@ -97,9 +97,16 @@ def __init__(self,
             self.decays.append(PyDarkNewsDecay(dec_case,
                                                table_dir = self.table_dir + table_subdirs))
 
-    def SaveCrossSectionTables(self):
+    def SaveCrossSectionTables(self,fill_tables_at_exit=True):
+        if not fill_tables_at_exit:
+            print('WARNING: Saving tables without filling PyDarkNewsCrossSection interpolation tables. Future updates to DarkNews can lead to inconsistent behavior if new entries are ever added to this table')
         for cross_section in self.cross_sections:
+            if fill_tables_at_exit:
+                print('Filling cross section table at %s'%cross_section.table_dir)
+                num = cross_section.FillInterpolationTables()
+                print('Added %d points'%num)
             cross_section.SaveInterpolationTables()
+
 
 
 
@@ -112,7 +119,7 @@ def __init__(self,
                  table_dir=None, # table to store 
                  tolerance=1e-6, # supposed to represent machine epsilon
                  interp_tolerance=5e-2, # relative interpolation tolerance
-                 interpolate_differential = False):
+                 interpolate_differential=False):
 
         DarkNewsCrossSection.__init__(self) # C++ constructor
 
@@ -170,23 +177,22 @@ def _redefine_interpolation_objects(self,total=False,diff=False):
                 # If we only have two energy points, don't try to construct interpolator
                 if len(np.unique(self.differential_cross_section_table[:,0])) <= 2: return
                 self.differential_cross_section_interpolator = CloughTocher2DInterpolator(self.differential_cross_section_table[:,:2],
-                                                                                        self.differential_cross_section_table[:,2],
-                                                                                        rescale=True)
+                                                                                          self.differential_cross_section_table[:,2],
+                                                                                          rescale=True)
 
     # Check whether we have close-enough entries in the intrepolation tables
-    def _query_interpolation_table(self,inputs,mode):
-        # 
-        # returns:
-        # 0 if we are not close enough to any points in the interpolation table
-        # otherwise, returns the desired interpolated value
-
+    def _interpolation_flags(self,inputs,mode):
+        #
+        # returns UseSinglePoint,Interpolate,closest_idx
+        # UseSinglePoint: whether to use a single point in table
+        # Interpolate: whether to interpolate bewteen different points
+        # closest_idx: index of closest point in table (for UseSinglePoint)
+
         # Determine which table we are using
         if mode=='total':
             interp_table = self.total_cross_section_table
-            interpolator = self.total_cross_section_interpolator
         elif mode=='differential':
             interp_table = self.differential_cross_section_table
-            interpolator = self.differential_cross_section_interpolator
         else:
             print('Invalid interpolation table mode %s'%mode)
             exit(0)
@@ -197,20 +203,15 @@ def _query_interpolation_table(self,inputs,mode):
         # bools to keep track of whether to use a single point or interpolate
         UseSinglePoint = True
         Interpolate = True
-        prev_closest_idx = None
+        # First check whether we have a close-enough single point
+        closest_idx = np.argmin(np.sum(np.abs(interp_table[:,:-1] - inputs)))
+        diff = (interp_table[closest_idx,:-1] - inputs)/inputs
+        if np.all(np.abs(diff)<self.tolerance): 
+            UseSinglePoint = True
         for i,input in enumerate(inputs):
             closest_idx = np.argmin(np.abs(interp_table[:,i] - input))
             diff = (interp_table[closest_idx,i] - input)/input # relative difference
-            # First check if we are close enough to use a single point
-            if np.abs(diff) >= self.tolerance:
-                # We are not close enough to use one existing table point
-                UseSinglePoint = False
-            elif UseSinglePoint:
-                # Check whether the closest_idx found here is the same as the previous
-                if i>0 and closest_idx != prev_closest_idx:
-                    UseSinglePoint = False
-                prev_closest_idx = closest_idx
-            # Now check if we are close enough to interpolate
+            # Check if we are close enough to interpolate
             if np.abs(diff) >= self.interp_tolerance:
                 Interpolate = False
             else:
@@ -239,13 +240,78 @@ def _query_interpolation_table(self,inputs,mode):
                         # check if the node above is also within the interpolation tolerance
                         if diff_above>=self.interp_tolerance:
                             Interpolate = False
+        return UseSinglePoint, Interpolate, closest_idx
+
+    # return entries in interpolation table if we have inputs 
+    def _query_interpolation_table(self,inputs,mode):
+        # 
+        # returns:
+        # 0 if we are not close enough to any points in the interpolation table
+        # otherwise, returns the desired interpolated value
+
+        # Determine which table we are using
+        if mode=='total':
+            interp_table = self.total_cross_section_table
+            interpolator = self.total_cross_section_interpolator
+        elif mode=='differential':
+            interp_table = self.differential_cross_section_table
+            interpolator = self.differential_cross_section_interpolator
+        else:
+            print('Invalid interpolation table mode %s'%mode)
+            exit(0)
+
+        UseSinglePoint, Interpolate, closest_idx = self._interpolation_flags(inputs,mode)
+
         if UseSinglePoint:
             return interp_table[closest_idx,-1]
         elif Interpolate:
             return interpolator(inputs)
         else:
             return 0    
 
+    # Fills the total and differential cross section tables within interp_tolerance
+    def FillInterpolationTables(self,total=True,diff=True):
+        Emin = (1.0+self.tolerance)*self.ups_case.Ethreshold
+        Emax = np.max(self.total_cross_section_table[:,0])
+        num_added_points = 0
+        if total:
+            E = Emin
+            while E<Emax:
+                UseSinglePoint,Interpolate,_ = self._interpolation_flags([E],'total')
+                if not (UseSinglePoint or Interpolate):
+                    xsec = self.ups_case.scalar_total_xsec(E)
+                    self.total_cross_section_table = np.append(self.total_cross_section_table,[[E,xsec]],axis=0)
+                    num_added_points+=1
+                E *= (1+self.interp_tolerance)
+        if diff:
+            # interaction record to calculate Q2 bounds
+            interaction = LI.dataclasses.InteractionRecord()
+            interaction.signature.primary_type = self.GetPossiblePrimaries()[0] # only one primary
+            interaction.signature.target_type = self.GetPossibleTargets()[0] # only one target
+            interaction.target_mass = self.ups_case.MA
+            E = Emin
+            zmin,zmax = self.tolerance,1
+            z = zmin
+            while E < Emax:
+                interaction.primary_momentum = [E,0,0,0]
+                Q2min = self.Q2Min(interaction)
+                Q2max = self.Q2Max(interaction)
+                z = zmin
+                while z < zmax:
+                    Q2 = Q2min + z*(Q2max-Q2min)
+                    UseSinglePoint,Interpolate,_ = self._interpolation_flags([E,z],'differential')
+                    if not (UseSinglePoint or Interpolate):
+                        dxsec = self.ups_case.diff_xsec_Q2(E, Q2).item()
+                        self.differential_cross_section_table = np.append(self.differential_cross_section_table,[[E,z,dxsec]],axis=0)
+                        num_added_points+=1
+                    z *= (1+self.interp_tolerance)
+                E *= (1+self.interp_tolerance)
+        self._redefine_interpolation_objects(total=total,diff=diff)
+        return num_added_points
+
+
+
+
     # Saves the tables for the scipy interpolation objects
     def SaveInterpolationTables(self,total=True,diff=True):
         if total:

resources/Examples/Example1/DIS_IceCube.py

@@ -0,0 +1,66 @@
+import numpy as np
+import os
+
+import leptoninjector as LI
+import sys
+sys.path.insert(1,'/n/holylfs05/LABS/arguelles_delgado_lab/Everyone/nkamp/LIV2/sources/LeptonInjector/python')
+from LIController import LIController
+
+LI_SRC = os.environ.get('LEPTONINJECTOR_SRC')
+
+# Number of events to inject
+events_to_inject = 1000
+
+# Expeirment to run
+experiment = 'IceCube'
+
+# Define the controller
+controller = LIController(events_to_inject,
+                          experiment)
+
+# Particle to inject
+primary_type = LI.dataclasses.Particle.ParticleType.NuMu
+
+# Cross Section Model
+xsfiledir = LI_SRC+'resources/CrossSectionTables/DISSplines/'
+target_type = LI.dataclasses.Particle.ParticleType.Nucleon
+
+DIS_xs = LI.crosssections.DISFromSpline(xsfiledir+'test_xs.fits',
+                                        xsfiledir+'test_xs_total.fits',
+                                        [primary_type],
+                                        [target_type])
+
+primary_xs = LI.crosssections.CrossSectionCollection(primary_type, [DIS_xs])
+controller.SetCrossSections(primary_xs)
+
+# Primary distributions
+primary_injection_distributions = {}
+primary_physical_distributions = {}
+
+# energy distribution
+edist = LI.distributions.PowerLaw(2,1e3,1e6)
+primary_injection_distributions['energy'] = edist
+primary_physical_distributions['energy'] = edist
+
+# direction distribution
+direction_distribution = LI.distributions.IsotropicDirection()
+primary_injection_distributions['direction'] = direction_distribution
+primary_physical_distributions['direction'] = direction_distribution
+
+# position distribution
+muon_range_func = LI.distributions.LeptonDepthFunction()
+position_distribution = LI.distributions.ColumnDepthPositionDistribution(600, 600.,
+                                                                         muon_range_func,
+                                                                         set(controller.GetEarthModelTargets()[0]))
+primary_injection_distributions['position'] = position_distribution
+
+# SetProcesses
+controller.SetProcesses(primary_type,
+                        primary_injection_distributions,
+                        primary_physical_distributions)
+
+controller.Initialize()
+
+events = controller.GenerateEvents()
+
+controller.SaveEvents('output/MINERvA_Dipole_M%2.2f_mu%2.2e_example.hdf5'%(model_kwargs['m4'],model_kwargs['mu_tr_mu4']))

resources/Examples/Example2/DipolePortal_MINERvA.py

@@ -0,0 +1,82 @@
+import numpy as np
+import os
+
+import leptoninjector as LI
+import sys
+sys.path.insert(1,'/n/holylfs05/LABS/arguelles_delgado_lab/Everyone/nkamp/LIV2/sources/LeptonInjector/python')
+from LIController import LIController
+
+LI_SRC = os.environ.get('LEPTONINJECTOR_SRC')
+
+# Define a DarkNews model 
+model_kwargs = {
+    'm4': 0.47,#0.140,
+    'mu_tr_mu4': 2.5e-6, #1e-6, # GeV^-1
+    'UD4': 0,
+    'Umu4': 0,
+    'epsilon': 0.0,
+    'gD': 0.0,
+    'decay_product':'photon',
+    'noHC':True,
+    'HNLtype':"dirac",
+}
+
+# Number of events to inject
+events_to_inject = 1000
+
+# Expeirment to run
+experiment = 'MINERvA'
+
+# Define the controller
+controller = LIController(events_to_inject,
+                          experiment)
+
+# Particle to inject
+primary_type = LI.dataclasses.Particle.ParticleType.NuMu
+
+# Define DarkNews Model
+table_dir = LI_SRC+'/resources/CrossSectionTables/DarkNewsTables/Dipole_M%2.2f_mu%2.2e/'%(model_kwargs['m4'],model_kwargs['mu_tr_mu4'])
+controller.InputDarkNewsModel(primary_type,
+                              table_dir,
+                              model_kwargs)
+
+# Primary distributions
+primary_injection_distributions = {}
+primary_physical_distributions = {}
+
+# energy distribution
+flux_file = LI_SRC + '/resources/Fluxes/NUMI_Flux_Tables/ME_FHC_numu.txt'
+edist = LI.distributions.TabulatedFluxDistribution(flux_file, True)
+edist_gen = LI.distributions.TabulatedFluxDistribution(1.05*model_kwargs['m4'], 10, flux_file, False)
+primary_injection_distributions['energy'] = edist_gen
+primary_physical_distributions['energy'] = edist
+
+# Flux normalization: go from cm^-2 to m^-2
+flux_units = LI.distributions.NormalizationConstant(1e4)
+primary_physical_distributions['normalization'] = flux_units
+
+# direction distribution
+direction_distribution = LI.distributions.FixedDirection(LI.math.Vector3D(0, 0, 1.0))
+primary_injection_distributions['direction'] = direction_distribution
+primary_physical_distributions['direction'] = direction_distribution
+
+# position distribution
+decay_range_func = LI.distributions.DecayRangeFunction(model_kwargs['m4'],
+                                                       controller.DN_min_decay_width,
+                                                       3,
+                                                       240)
+position_distribution = LI.distributions.RangePositionDistribution(1.24, 5.,
+                                                                   decay_range_func,
+                                                                   set(controller.GetEarthModelTargets()[0]))
+primary_injection_distributions['position'] = position_distribution
+
+# SetProcesses
+controller.SetProcesses(primary_type,
+                        primary_injection_distributions,
+                        primary_physical_distributions)
+
+controller.Initialize()
+
+events = controller.GenerateEvents()
+
+controller.SaveEvents('output/MINERvA_Dipole_M%2.2f_mu%2.2e_example.hdf5'%(model_kwargs['m4'],model_kwargs['mu_tr_mu4']))

resources/Examples/Example2/DipolePortal_MiniBooNE.py

@@ -23,7 +23,7 @@
 }
 
 # Number of events to inject
-events_to_inject = 10
+events_to_inject = 1000
 
 # Expeirment to run
 experiment = 'MiniBooNE'
@@ -80,38 +80,5 @@
 
 events = controller.GenerateEvents()
 
-controller.SaveEvents('Dipole_M%2.2f_mu%2.2e_example.hdf5'%(model_kwargs['m4'],model_kwargs['mu_tr_mu4']))
-
-
-# Analysis
-
-upscattering_vertices = {}
-targets = []
-decay_vertices = []
-for tree in events:
-    for datum in tree.tree:
-        if datum.record.signature.primary_type == LI.dataclasses.Particle.ParticleType.NuMu:
-            target = datum.record.signature.target_type
-            if target not in upscattering_vertices.keys():
-                upscattering_vertices[target] = []
-            upscattering_vertices[target].append(datum.record.interaction_vertex)
-        else:
-            decay_vertices.append(datum.record.interaction_vertex)
-
-for k,vertices in upscattering_vertices.items():
-    vertices = np.array(vertices)
-    plt.scatter(vertices[:,2],vertices[:,0],label=k,alpha=0.3)
-plt.xlabel('Z [m]')
-plt.ylabel('X [m]')
-plt.legend()
-plt.show()
-plt.savefig('Figures/%s_upscattering.png'%experiment,dpi=100)
-
-decay_vertices = np.array(decay_vertices)
-plt.scatter(decay_vertices[:,2],decay_vertices[:,0],alpha=0.3)
-plt.xlabel('Z [m]')
-plt.ylabel('X [m]')
-plt.legend()
-plt.show()
-plt.savefig('Figures/%s_decay.png'%experiment,dpi=100)
+controller.SaveEvents('output/MiniBooNE_Dipole_M%2.2f_mu%2.2e_example.hdf5'%(model_kwargs['m4'],model_kwargs['mu_tr_mu4']))