Module xenodiffusionscope.LCEPattern

Expand source code
import matplotlib.pyplot as plt
import numpy as np
import scipy.integrate
import scipy.interpolate
from matplotlib.patches import Circle

from .TPC import TPC

#Matplotlib plots should be handled in separate file with Plotter.

class LCEPattern:
    '''
    Ad-hoc LCE maps. From pos and number of electrons to photons detected.
    '''
    
    def __init__(self, TPC):
        
        # TO BE CLEANED UP
        self.TPC = TPC
        self.z_anode = TPC.z_anode
        self.r_max = TPC.radius
        # Since some sensors stick out of the array, we need to consider
        # photons ending outside the nominal radius of the TPC in the pattern.
        self.r_pattern_max = 100
        
        #self.load_patterns()
    
    def define_pattern_props(self, x_bin_step, y_bin_step, n_traces, 
                             smooth_pattern, force_traces = False):
        '''
        Call this function to define the pattern calculation properties.
        '''
        
        self.x_bin_step = x_bin_step #mm #to put in config
        self.y_bin_step = y_bin_step #mm #to put in config
        self.pattern_bin_area = self.x_bin_step * self.y_bin_step
        if force_traces != True:
            assert n_traces < 1e7, 'Asked for too many traces. Proceed on your own accord with force_traces = True'
        self.n_traces = int(n_traces)
        self.smooth_pattern = smooth_pattern
        return None
    
    @classmethod
    def get_xy_on_plane(cls,x0,y0,z0,directions,z_prime):
        '''
        Return the (x,y) of the intersection of a vector with direction 
        (theta,phi), physics convention, starting at point (x0,y0,z0) and
        finishing at the plane z=z_prime.
        '''
        x_prime = x0 + (z_prime-z0) * np.cos(directions[:,1]) * np.tan(directions[:,0])
        y_prime = y0 + (z_prime-z0) * np.sin(directions[:,1]) * np.tan(directions[:,0])
        final_pos = np.stack((x_prime, y_prime), axis = 1)
        return final_pos
    
    def get_xy_on_circ_array(self, final_pos):
        '''Return the x,y positions of the electron cluds ending inside the 
        radius to considere for the pattern (TPC + sensors outside).'''
        final_r = TPC.get_r(final_pos[:,0],final_pos[:,1])
        mask_r = final_r < self.r_pattern_max
        final_pos_array = final_pos[mask_r]
        return final_pos_array
    
    def get_hits_on_circ_array(self, x0,y0,z0):
        '''
        Get the positions of the toys that hit the circular area of the array.
        '''
        
        thetas = np.random.uniform(0,np.pi/2,self.n_traces)
        phis = np.random.uniform(0,2*np.pi,self.n_traces)
        directions = np.stack((thetas,phis), axis = 1)
        
        final_pos = self.get_xy_on_plane(x0,y0,z0,directions, self.z_anode)
        final_pos_array = self.get_xy_on_circ_array(final_pos)
        return final_pos_array
    

    def print_stats_of_hits(self, hits, n_traces):
        '''Print statistics of the hits.'''

        print('Initial number of photons: %s'%n_traces)
        print('Number of photons hit: '+
              '%s (%.2f%% of produced, %.2f of '+
              'full emission)'%(len(hits),
                                len(hits)/n_traces*100,
                                len(hits)/n_traces*100/2))
        return None
    
    def get_pattern_density_hist2d(self, pos):
        '''Calculate the density of photons in x,y bins. Result is a 2D 
        histogram, already normalised to the bin area and number of initial 
        photons.'''
        
        x_min = y_min = -np.ceil(self.r_pattern_max*1.1)
        x_max = y_max = np.ceil(self.r_pattern_max*1.1)
        
        x_bin_sides = np.arange(x_min,x_max+self.x_bin_step*0.1,self.x_bin_step)
        y_bin_sides = np.arange(y_min,y_max+self.x_bin_step*0.1,self.y_bin_step)
        
        x_bin_middles = x_bin_sides[:-1] + self.x_bin_step/2
        y_bin_middles = y_bin_sides[:-1] + self.y_bin_step/2
        
        hist2d = np.histogram2d(pos[:,0],pos[:,1],
                        bins = (x_bin_sides,y_bin_sides),
                                density = False)
        #density could be True but we can also do it by hand given the bin area 
        #and total of photons. It's even better because it can be normalized properly 
        #from the start taking into account the photonos projected downwards and the
        #ones that miss the array.
        assert np.sum(hist2d[0]) == len(pos), ('Lost some photons on the histogram??\n'+
                                               'N_toys: %d\n' %len(pos)+
                                               'N_hist2d: %d' %len(hist2d[2]))
        
        
        #total_in_hist = np.sum(hist2d[0])
        #double the traces innormalization because of under side of sphere
        hist_density = hist2d[0]/self.pattern_bin_area/(self.n_traces*2) #fraction/mm^2;
        return x_bin_middles,y_bin_middles,hist_density
    
    def make_pattern_density(self,pos):
        '''
        Takes the toy-MC results and makes the 2D density histogram, 
        normalized to the total number of traces produced and bin area.
        Returns either the interpolative function.
        Params:
          * pos: (N,2) array of the positions of hits in the circular array
          * x_bin_step: size of the hist2d bin length in x
          * y_bin_step: size of the hist2d bin length in y
        '''
        if self.smooth_pattern:
            s = 0.0000001
        else:
            s = 0
        x_bin_middles,y_bin_middles,hist_density = self.get_pattern_density_hist2d(pos)
        interp2s = scipy.interpolate.RectBivariateSpline(x_bin_middles, 
                                                         y_bin_middles,
                                                         hist_density,
                                                         s = s)
        return interp2s
    
    def make_pattern_from_pos(self,x0,y0,z0):
        '''Merger function to get the pattern from a given x,y,z position.'''
        toy_events = self.get_hits_on_circ_array(x0,y0,z0)
        pattern = self.make_pattern_density(toy_events)
        return pattern

Classes

class LCEPattern (TPC)

Ad-hoc LCE maps. From pos and number of electrons to photons detected.

Expand source code
class LCEPattern:
    '''
    Ad-hoc LCE maps. From pos and number of electrons to photons detected.
    '''
    
    def __init__(self, TPC):
        
        # TO BE CLEANED UP
        self.TPC = TPC
        self.z_anode = TPC.z_anode
        self.r_max = TPC.radius
        # Since some sensors stick out of the array, we need to consider
        # photons ending outside the nominal radius of the TPC in the pattern.
        self.r_pattern_max = 100
        
        #self.load_patterns()
    
    def define_pattern_props(self, x_bin_step, y_bin_step, n_traces, 
                             smooth_pattern, force_traces = False):
        '''
        Call this function to define the pattern calculation properties.
        '''
        
        self.x_bin_step = x_bin_step #mm #to put in config
        self.y_bin_step = y_bin_step #mm #to put in config
        self.pattern_bin_area = self.x_bin_step * self.y_bin_step
        if force_traces != True:
            assert n_traces < 1e7, 'Asked for too many traces. Proceed on your own accord with force_traces = True'
        self.n_traces = int(n_traces)
        self.smooth_pattern = smooth_pattern
        return None
    
    @classmethod
    def get_xy_on_plane(cls,x0,y0,z0,directions,z_prime):
        '''
        Return the (x,y) of the intersection of a vector with direction 
        (theta,phi), physics convention, starting at point (x0,y0,z0) and
        finishing at the plane z=z_prime.
        '''
        x_prime = x0 + (z_prime-z0) * np.cos(directions[:,1]) * np.tan(directions[:,0])
        y_prime = y0 + (z_prime-z0) * np.sin(directions[:,1]) * np.tan(directions[:,0])
        final_pos = np.stack((x_prime, y_prime), axis = 1)
        return final_pos
    
    def get_xy_on_circ_array(self, final_pos):
        '''Return the x,y positions of the electron cluds ending inside the 
        radius to considere for the pattern (TPC + sensors outside).'''
        final_r = TPC.get_r(final_pos[:,0],final_pos[:,1])
        mask_r = final_r < self.r_pattern_max
        final_pos_array = final_pos[mask_r]
        return final_pos_array
    
    def get_hits_on_circ_array(self, x0,y0,z0):
        '''
        Get the positions of the toys that hit the circular area of the array.
        '''
        
        thetas = np.random.uniform(0,np.pi/2,self.n_traces)
        phis = np.random.uniform(0,2*np.pi,self.n_traces)
        directions = np.stack((thetas,phis), axis = 1)
        
        final_pos = self.get_xy_on_plane(x0,y0,z0,directions, self.z_anode)
        final_pos_array = self.get_xy_on_circ_array(final_pos)
        return final_pos_array
    

    def print_stats_of_hits(self, hits, n_traces):
        '''Print statistics of the hits.'''

        print('Initial number of photons: %s'%n_traces)
        print('Number of photons hit: '+
              '%s (%.2f%% of produced, %.2f of '+
              'full emission)'%(len(hits),
                                len(hits)/n_traces*100,
                                len(hits)/n_traces*100/2))
        return None
    
    def get_pattern_density_hist2d(self, pos):
        '''Calculate the density of photons in x,y bins. Result is a 2D 
        histogram, already normalised to the bin area and number of initial 
        photons.'''
        
        x_min = y_min = -np.ceil(self.r_pattern_max*1.1)
        x_max = y_max = np.ceil(self.r_pattern_max*1.1)
        
        x_bin_sides = np.arange(x_min,x_max+self.x_bin_step*0.1,self.x_bin_step)
        y_bin_sides = np.arange(y_min,y_max+self.x_bin_step*0.1,self.y_bin_step)
        
        x_bin_middles = x_bin_sides[:-1] + self.x_bin_step/2
        y_bin_middles = y_bin_sides[:-1] + self.y_bin_step/2
        
        hist2d = np.histogram2d(pos[:,0],pos[:,1],
                        bins = (x_bin_sides,y_bin_sides),
                                density = False)
        #density could be True but we can also do it by hand given the bin area 
        #and total of photons. It's even better because it can be normalized properly 
        #from the start taking into account the photonos projected downwards and the
        #ones that miss the array.
        assert np.sum(hist2d[0]) == len(pos), ('Lost some photons on the histogram??\n'+
                                               'N_toys: %d\n' %len(pos)+
                                               'N_hist2d: %d' %len(hist2d[2]))
        
        
        #total_in_hist = np.sum(hist2d[0])
        #double the traces innormalization because of under side of sphere
        hist_density = hist2d[0]/self.pattern_bin_area/(self.n_traces*2) #fraction/mm^2;
        return x_bin_middles,y_bin_middles,hist_density
    
    def make_pattern_density(self,pos):
        '''
        Takes the toy-MC results and makes the 2D density histogram, 
        normalized to the total number of traces produced and bin area.
        Returns either the interpolative function.
        Params:
          * pos: (N,2) array of the positions of hits in the circular array
          * x_bin_step: size of the hist2d bin length in x
          * y_bin_step: size of the hist2d bin length in y
        '''
        if self.smooth_pattern:
            s = 0.0000001
        else:
            s = 0
        x_bin_middles,y_bin_middles,hist_density = self.get_pattern_density_hist2d(pos)
        interp2s = scipy.interpolate.RectBivariateSpline(x_bin_middles, 
                                                         y_bin_middles,
                                                         hist_density,
                                                         s = s)
        return interp2s
    
    def make_pattern_from_pos(self,x0,y0,z0):
        '''Merger function to get the pattern from a given x,y,z position.'''
        toy_events = self.get_hits_on_circ_array(x0,y0,z0)
        pattern = self.make_pattern_density(toy_events)
        return pattern

Static methods

def get_xy_on_plane(x0, y0, z0, directions, z_prime)

Return the (x,y) of the intersection of a vector with direction (theta,phi), physics convention, starting at point (x0,y0,z0) and finishing at the plane z=z_prime.

Expand source code
@classmethod
def get_xy_on_plane(cls,x0,y0,z0,directions,z_prime):
    '''
    Return the (x,y) of the intersection of a vector with direction 
    (theta,phi), physics convention, starting at point (x0,y0,z0) and
    finishing at the plane z=z_prime.
    '''
    x_prime = x0 + (z_prime-z0) * np.cos(directions[:,1]) * np.tan(directions[:,0])
    y_prime = y0 + (z_prime-z0) * np.sin(directions[:,1]) * np.tan(directions[:,0])
    final_pos = np.stack((x_prime, y_prime), axis = 1)
    return final_pos

Methods

def define_pattern_props(self, x_bin_step, y_bin_step, n_traces, smooth_pattern, force_traces=False)

Call this function to define the pattern calculation properties.

Expand source code
def define_pattern_props(self, x_bin_step, y_bin_step, n_traces, 
                         smooth_pattern, force_traces = False):
    '''
    Call this function to define the pattern calculation properties.
    '''
    
    self.x_bin_step = x_bin_step #mm #to put in config
    self.y_bin_step = y_bin_step #mm #to put in config
    self.pattern_bin_area = self.x_bin_step * self.y_bin_step
    if force_traces != True:
        assert n_traces < 1e7, 'Asked for too many traces. Proceed on your own accord with force_traces = True'
    self.n_traces = int(n_traces)
    self.smooth_pattern = smooth_pattern
    return None
def get_hits_on_circ_array(self, x0, y0, z0)

Get the positions of the toys that hit the circular area of the array.

Expand source code
def get_hits_on_circ_array(self, x0,y0,z0):
    '''
    Get the positions of the toys that hit the circular area of the array.
    '''
    
    thetas = np.random.uniform(0,np.pi/2,self.n_traces)
    phis = np.random.uniform(0,2*np.pi,self.n_traces)
    directions = np.stack((thetas,phis), axis = 1)
    
    final_pos = self.get_xy_on_plane(x0,y0,z0,directions, self.z_anode)
    final_pos_array = self.get_xy_on_circ_array(final_pos)
    return final_pos_array
def get_pattern_density_hist2d(self, pos)

Calculate the density of photons in x,y bins. Result is a 2D histogram, already normalised to the bin area and number of initial photons.

Expand source code
def get_pattern_density_hist2d(self, pos):
    '''Calculate the density of photons in x,y bins. Result is a 2D 
    histogram, already normalised to the bin area and number of initial 
    photons.'''
    
    x_min = y_min = -np.ceil(self.r_pattern_max*1.1)
    x_max = y_max = np.ceil(self.r_pattern_max*1.1)
    
    x_bin_sides = np.arange(x_min,x_max+self.x_bin_step*0.1,self.x_bin_step)
    y_bin_sides = np.arange(y_min,y_max+self.x_bin_step*0.1,self.y_bin_step)
    
    x_bin_middles = x_bin_sides[:-1] + self.x_bin_step/2
    y_bin_middles = y_bin_sides[:-1] + self.y_bin_step/2
    
    hist2d = np.histogram2d(pos[:,0],pos[:,1],
                    bins = (x_bin_sides,y_bin_sides),
                            density = False)
    #density could be True but we can also do it by hand given the bin area 
    #and total of photons. It's even better because it can be normalized properly 
    #from the start taking into account the photonos projected downwards and the
    #ones that miss the array.
    assert np.sum(hist2d[0]) == len(pos), ('Lost some photons on the histogram??\n'+
                                           'N_toys: %d\n' %len(pos)+
                                           'N_hist2d: %d' %len(hist2d[2]))
    
    
    #total_in_hist = np.sum(hist2d[0])
    #double the traces innormalization because of under side of sphere
    hist_density = hist2d[0]/self.pattern_bin_area/(self.n_traces*2) #fraction/mm^2;
    return x_bin_middles,y_bin_middles,hist_density
def get_xy_on_circ_array(self, final_pos)

Return the x,y positions of the electron cluds ending inside the radius to considere for the pattern (TPC + sensors outside).

Expand source code
def get_xy_on_circ_array(self, final_pos):
    '''Return the x,y positions of the electron cluds ending inside the 
    radius to considere for the pattern (TPC + sensors outside).'''
    final_r = TPC.get_r(final_pos[:,0],final_pos[:,1])
    mask_r = final_r < self.r_pattern_max
    final_pos_array = final_pos[mask_r]
    return final_pos_array
def make_pattern_density(self, pos)

Takes the toy-MC results and makes the 2D density histogram, normalized to the total number of traces produced and bin area. Returns either the interpolative function.

Params

  • pos: (N,2) array of the positions of hits in the circular array
  • x_bin_step: size of the hist2d bin length in x
  • y_bin_step: size of the hist2d bin length in y
Expand source code
def make_pattern_density(self,pos):
    '''
    Takes the toy-MC results and makes the 2D density histogram, 
    normalized to the total number of traces produced and bin area.
    Returns either the interpolative function.
    Params:
      * pos: (N,2) array of the positions of hits in the circular array
      * x_bin_step: size of the hist2d bin length in x
      * y_bin_step: size of the hist2d bin length in y
    '''
    if self.smooth_pattern:
        s = 0.0000001
    else:
        s = 0
    x_bin_middles,y_bin_middles,hist_density = self.get_pattern_density_hist2d(pos)
    interp2s = scipy.interpolate.RectBivariateSpline(x_bin_middles, 
                                                     y_bin_middles,
                                                     hist_density,
                                                     s = s)
    return interp2s
def make_pattern_from_pos(self, x0, y0, z0)

Merger function to get the pattern from a given x,y,z position.

Expand source code
def make_pattern_from_pos(self,x0,y0,z0):
    '''Merger function to get the pattern from a given x,y,z position.'''
    toy_events = self.get_hits_on_circ_array(x0,y0,z0)
    pattern = self.make_pattern_density(toy_events)
    return pattern
def print_stats_of_hits(self, hits, n_traces)

Print statistics of the hits.

Expand source code
def print_stats_of_hits(self, hits, n_traces):
    '''Print statistics of the hits.'''

    print('Initial number of photons: %s'%n_traces)
    print('Number of photons hit: '+
          '%s (%.2f%% of produced, %.2f of '+
          'full emission)'%(len(hits),
                            len(hits)/n_traces*100,
                            len(hits)/n_traces*100/2))
    return None