SndlhcMuonReco.hough Class Reference

Public Member Functions
	__init__ (self, n_yH, yH_range, n_xH, xH_range, z_offset, Hformat, space_scale, det_Zlen, squaretheta=False, smooth=True)
	fit (self, hit_collection, is_scaled, draw, weights=None)
	fit_randomize (self, hit_collection, hit_d, n_random, is_scaled, draw, weights=None)

Public Attributes
	n_yH
	n_xH
	yH_range
	xH_range
	z_offset
	HoughSpace_format
	space_scale
	det_Zlen
	smooth
	yH_bins
	xH_bins
	cos_thetas
	sin_thetas
	xH_i
	n_yH_scaled
	n_xH_scaled
	yH_bins_scaled
	xH_bins_scaled
	cos_thetas_scaled
	sin_thetas_scaled
	xH_i_scaled
	accumulator

Detailed Description

 Hough transform implementation 

Definition at line 25 of file SndlhcMuonReco.py.

Constructor & Destructor Documentation

__init__()

SndlhcMuonReco.hough.__init__	(		self,
			n_yH,
			yH_range,
			n_xH,
			xH_range,
			z_offset,
			Hformat,
			space_scale,
			det_Zlen,
			squaretheta = False,
			smooth = True
	)

Definition at line 28 of file SndlhcMuonReco.py.

    def __init__(self, n_yH, yH_range, n_xH, xH_range, z_offset, Hformat, space_scale, det_Zlen, squaretheta = False, smooth = True) :
 
        self.n_yH = n_yH
        self.n_xH = n_xH
 
        self.yH_range = yH_range
        self.xH_range = xH_range
 
        self.z_offset = z_offset
        self.HoughSpace_format = Hformat
        self.space_scale = space_scale
        
        self.det_Zlen = det_Zlen
 
        self.smooth = smooth
 
        self.yH_bins = np.linspace(self.yH_range[0], self.yH_range[1], n_yH)
        if not squaretheta :
            self.xH_bins = np.linspace(self.xH_range[0], self.xH_range[1], n_xH)
        else :
            self.xH_bins = np.linspace(np.sign(self.xH_range[0])*(self.xH_range[0]**0.5), np.sign(self.xH_range[1])*(self.xH_range[1]**0.5), n_xH)
            self.xH_bins = np.sign(self.xH_bins)*np.square(self.xH_bins)
 
        self.cos_thetas = np.cos(self.xH_bins)
        self.sin_thetas = np.sin(self.xH_bins)
        
        self.xH_i = np.array(list(range(n_xH)))
 
        # A back-up Hough space designed to have more/less bins wrt the default one above.
        # It is useful when fitting some low-E muon tracks, which are curved due to mult. scattering.
        self.n_yH_scaled = int(n_yH*space_scale)
        self.n_xH_scaled = int(n_xH*space_scale)
        self.yH_bins_scaled = np.linspace(self.yH_range[0], self.yH_range[1], self.n_yH_scaled)
        if not squaretheta :
            self.xH_bins_scaled= np.linspace(self.xH_range[0], self.xH_range[1], self.n_xH_scaled)
        else :
            self.xH_bins_scaled = np.linspace(np.sign(self.xH_range[0])*(self.xH_range[0]**0.5), np.sign(self.xH_range[1])*(self.xH_range[1]**0.5), self.n_xH_scaled)
            self.xH_bins_scaled = np.sign(self.xH_bins_scaled)*np.square(self.xH_bins_scaled)
 
        self.cos_thetas_scaled = np.cos(self.xH_bins_scaled)
        self.sin_thetas_scaled = np.sin(self.xH_bins_scaled)
 
        self.xH_i_scaled = np.array(list(range(self.n_xH_scaled)))
 

Member Function Documentation

fit()

SndlhcMuonReco.hough.fit	(		self,
			hit_collection,
			is_scaled,
			draw,
			weights = None
	)

Definition at line 72 of file SndlhcMuonReco.py.

    def fit(self, hit_collection, is_scaled, draw, weights = None) :
 
        if not is_scaled:
           n_xH = self.n_xH
           n_yH = self.n_yH
           cos_thetas = self.cos_thetas
           sin_thetas = self.sin_thetas
           xH_bins = self.xH_bins
           yH_bins = self.yH_bins
           xH_i = self.xH_i
           res = self.res
        else:
           n_xH = self.n_xH_scaled
           n_yH = self.n_yH_scaled
           cos_thetas = self.cos_thetas_scaled
           sin_thetas = self.sin_thetas_scaled
           xH_bins = self.xH_bins_scaled
           yH_bins = self.yH_bins_scaled
           xH_i = self.xH_i_scaled
           res = self.res*self.space_scale
 
        self.accumulator = np.zeros((n_yH, n_xH))
        for i_hit, hit in enumerate(hit_collection) :
            shifted_hitZ = hit[0] - self.z_offset
            if self.HoughSpace_format == 'normal':
                 hit_yH = shifted_hitZ*cos_thetas + hit[1]*sin_thetas
            elif self.HoughSpace_format == 'linearSlopeIntercept':
                 hit_yH = hit[1] - shifted_hitZ*xH_bins
            elif self.HoughSpace_format== 'linearIntercepts':
                 hit_yH = (self.det_Zlen*hit[1] - shifted_hitZ*xH_bins)/(self.det_Zlen - shifted_hitZ)
            out_of_range = np.logical_and(hit_yH > self.yH_range[0], hit_yH < self.yH_range[1]) 
            hit_yH_i = np.floor((hit_yH[out_of_range] - self.yH_range[0])/(self.yH_range[1] - self.yH_range[0])*n_yH).astype(np.int_)
 
            if weights is not None :
                self.accumulator[hit_yH_i, xH_i[out_of_range]] += weights[i_hit]
            else :
                self.accumulator[hit_yH_i, xH_i[out_of_range]] += 1
 
        # Smooth accumulator
        if self.smooth_full :
            self.accumulator = scipy.ndimage.gaussian_filter(self.accumulator, self.sigma, truncate=self.truncate)
 
        # This might be useful for debugging, but leave out for now.
        if draw :
            plt.figure()
            plt.imshow(self.accumulator, origin = "lower", extent = [self.xH_range[0], self.xH_range[-1], self.yH_range[0], self.yH_range[-1]], aspect = "auto")
            if self.HoughSpace_format == 'normal':
               plt.xlabel(r"$\theta$ [rad]")
               plt.ylabel("r [cm]")
            elif self.HoughSpace_format == 'linearSlopeIntercept':
               plt.xlabel("slope")
               plt.ylabel("intercept @ 1st plane [cm]")
            elif self.HoughSpace_format == 'linearIntercepts':
               plt.xlabel("intercept @ last plane [cm]")
               plt.ylabel("intercept @ 1st plane [cm]")
            plt.tight_layout()
            plt.show()
 
        if self.smooth_full:
          # In case of multiple occurrences of the maximum values, argmax returns
          # the indices corresponding to the first occurrence(along 1st axis).
          i_max = np.unravel_index(self.accumulator.argmax(), self.accumulator.shape)
        else:
          # In case there are more than 1 bins with the maximal Nentries, check if the n-th quantile of
          # found peaks along 'slope' axis(yH axis) enloses <n-th quantile> portion of all maxima 'slope' bins
          # within a specified range. It is advisable that that range is consistent with
          # the detector angular resolution. 
          maxima = np.argwhere(self.accumulator == np.amax(self.accumulator))
          if len(maxima) == 1:
            i_max = maxima[0]
          elif len(maxima)==0 or (len(maxima) > 1 and self.HoughSpace_format == 'linearIntercepts'):
               # if no reasonable way to select btw maxima, force track-build failure
               # to be decided what to do in 'linearIntercepts' case and multiple maxima
               return(-999, -999)
          elif len(maxima) > 1 and abs(min([x[1] for x in maxima]) - max([x[1] for x in maxima])) < res:
               i_max = maxima[0]
          else:
            # FIXME: the next two lines can be a single line command for sure
            maxima_slopesAxis_list = list([x[1] for x in maxima])
            maxima_slopesAxis_list = np.asarray(maxima_slopesAxis_list)
            quantile = np.quantile(maxima_slopesAxis_list, self.n_quantile)
            Nwithin = 0
            # FIXME: a more elegant 'hidden' loop is maybe possible here too
            for item in maxima:
               if abs(item[1]-quantile)< res:
                 Nwithin += 1
            if Nwithin/len(maxima) > self.n_quantile: 
               i_x = min([x[1] for x in maxima], key=lambda b: abs(b-quantile))
               for im in maxima:
                 if im[1] == i_x : i_max = im
            else: 
                 # if no reasonable way to select btw maxima, force track-build failure
                 return(-999, -999)
 
        found_yH = yH_bins[int(i_max[0])]
        found_xH = xH_bins[int(i_max[1])]
        
        if self.HoughSpace_format == 'normal':
           slope = -1./np.tan(found_xH)
           interceptShift = found_yH/np.sin(found_xH)
           intercept = (np.tan(found_xH)*interceptShift + self.z_offset)/np.tan(found_xH)
        elif self.HoughSpace_format == 'linearSlopeIntercept':
           slope = found_xH
           intercept = found_yH - slope*self.z_offset
        elif self.HoughSpace_format == 'linearIntercepts':
           slope = (found_xH - found_yH)/self.det_Zlen
           intercept = found_yH - slope*self.z_offset
        
        return (slope, intercept)
 

fit_randomize()

SndlhcMuonReco.hough.fit_randomize	(		self,
			hit_collection,
			hit_d,
			n_random,
			is_scaled,
			draw,
			weights = None
	)

Definition at line 182 of file SndlhcMuonReco.py.

    def fit_randomize(self, hit_collection, hit_d, n_random, is_scaled, draw, weights = None) :
        success = True
        if not len(hit_collection) :
            return (-1, -1, [[],[]], [], False)
 
        # Randomize hits
        if (n_random > 0) :
            random_hit_collection = []
            for i_random in range(n_random) :
                random_hits = np.random.uniform(size = hit_collection.shape) - 0.5
                random_hits *= hit_d
                random_hits += hit_collection
                random_hit_collection.append(random_hits)
 
            random_hit_collection = np.concatenate(random_hit_collection)
            if weights is not None :
                weights = np.tile(weights, n_random)
 
            fit = self.fit(random_hit_collection, is_scaled, draw, weights)
        else :
            fit = self.fit(hit_collection, is_scaled, draw, weights)
 
        return fit
 

Member Data Documentation

accumulator

SndlhcMuonReco.hough.accumulator

Definition at line 93 of file SndlhcMuonReco.py.

cos_thetas

SndlhcMuonReco.hough.cos_thetas

Definition at line 51 of file SndlhcMuonReco.py.

cos_thetas_scaled

SndlhcMuonReco.hough.cos_thetas_scaled

Definition at line 67 of file SndlhcMuonReco.py.

det_Zlen

SndlhcMuonReco.hough.det_Zlen

Definition at line 40 of file SndlhcMuonReco.py.

HoughSpace_format

SndlhcMuonReco.hough.HoughSpace_format

Definition at line 37 of file SndlhcMuonReco.py.

n_xH

SndlhcMuonReco.hough.n_xH

Definition at line 31 of file SndlhcMuonReco.py.

n_xH_scaled

SndlhcMuonReco.hough.n_xH_scaled

Definition at line 59 of file SndlhcMuonReco.py.

n_yH

SndlhcMuonReco.hough.n_yH

Definition at line 30 of file SndlhcMuonReco.py.

n_yH_scaled

SndlhcMuonReco.hough.n_yH_scaled

Definition at line 58 of file SndlhcMuonReco.py.

sin_thetas

SndlhcMuonReco.hough.sin_thetas

Definition at line 52 of file SndlhcMuonReco.py.

sin_thetas_scaled

SndlhcMuonReco.hough.sin_thetas_scaled

Definition at line 68 of file SndlhcMuonReco.py.

smooth

SndlhcMuonReco.hough.smooth

Definition at line 42 of file SndlhcMuonReco.py.

space_scale

SndlhcMuonReco.hough.space_scale

Definition at line 38 of file SndlhcMuonReco.py.

xH_bins

SndlhcMuonReco.hough.xH_bins

Definition at line 46 of file SndlhcMuonReco.py.

xH_bins_scaled

SndlhcMuonReco.hough.xH_bins_scaled

Definition at line 62 of file SndlhcMuonReco.py.

xH_i

SndlhcMuonReco.hough.xH_i

Definition at line 54 of file SndlhcMuonReco.py.

xH_i_scaled

SndlhcMuonReco.hough.xH_i_scaled

Definition at line 70 of file SndlhcMuonReco.py.

xH_range

SndlhcMuonReco.hough.xH_range

Definition at line 34 of file SndlhcMuonReco.py.

yH_bins

SndlhcMuonReco.hough.yH_bins

Definition at line 44 of file SndlhcMuonReco.py.

yH_bins_scaled

SndlhcMuonReco.hough.yH_bins_scaled

Definition at line 60 of file SndlhcMuonReco.py.

yH_range

SndlhcMuonReco.hough.yH_range

Definition at line 33 of file SndlhcMuonReco.py.

z_offset

SndlhcMuonReco.hough.z_offset

Definition at line 36 of file SndlhcMuonReco.py.

---

The documentation for this class was generated from the following file:

• python/SndlhcMuonReco.py