eic-opticks/ana/toymuon.py

0001 #!/usr/bin/env python
0002
0003 import numpy as np
0004 import math
0005 from scipy.spatial.transform import Rotation as R
0006
0007 m = R.from_euler('z', 90, degrees=True)
0008
0009 try:
0010     import pyvista as pv
0011 except ImportError:
0012     pv = None
0013 pass
0014
0015 try:
0016     import matplotlib.pyplot as plt
0017 except ImportError:
0018     plt = None
0019 pass
0020
0021
0022 def yan_0():
0023
0024     rndmcostheta = random.uniform(-1,1);
0025     rndmphi = random.uniform(0,math.pi*2);
0026     r = math.sqrt(1-rndmcostheta*rndmcostheta);
0027     x = r*math.cos(rndmphi);
0028     y = r*math.sin(rndmphi);
0029
0030     return x,y,z
0031
0032
0033 def uniform_points_on_sphere(size=1000):
0034
0035     u0 = np.random.uniform(low=-1, high=1, size=size)
0036     u1 = np.random.uniform(low=0, high=np.pi*2., size=size)
0037
0038     pos = np.zeros( (size,3) , dtype=np.float32)
0039
0040     azimuth = u1
0041     cos_polar = u0
0042     sin_polar = np.sqrt( 1.-u0*u0 )
0043
0044     x = sin_polar*np.cos(azimuth)
0045     y = sin_polar*np.sin(azimuth)
0046     z = cos_polar
0047
0048     pos[:,0] = x
0049     pos[:,1] = y
0050     pos[:,2] = z
0051
0052     return pos
0053
0054
0055 class Generator(object):
0056     def plot(self, ax):
0057         ax.plot(self.x,self.y)
0058         ax.plot(self.x,self.cdf)
0059
0060     def __call__(self, u ):
0061         """
0062         :param u: uniform random number in range 0->1 (or numpy array of such randoms)
0063         :return x: generated sample that follows the origin pdf by interpolation of inverse CDF
0064         """
0065         return np.interp( u, self.idom, self.icdf )
0066
0067
0068 class MuonZenithAngle(Generator):
0069     """
0070     Googling for "Gaisser formula plot" yields:
0071
0072     "Muon Simulation at the Daya Bay Site",  Mengyun, Guan (2010)
0073
0074     https://escholarship.org/uc/item/6jm8g76d
0075
0076     Fig 10 from the above paper looks like a bit like   cos^2( polar-pi/4 ) - 0.5
0077     so lets use that as an ad-hoc approximation.
0078
0079     For a review on "COSMIC RAY MUON PHYSICS" see
0080
0081     * https://core.ac.uk/download/pdf/25279944.pdf
0082     """
0083
0084     def __init__(self, x, num_idom=1000):
0085
0086         xp = x-0.25*np.pi
0087         y = np.cos(xp)*np.cos(xp) - 0.5
0088         # adhoc choice of function, should really use full Gaisser with fitted params from measuremnts
0089         #y /= y.sum()
0090
0091         cdf = np.cumsum(y)
0092         cdf /= cdf[-1]       # CDFs tend to 1.
0093
0094         dom = x                             # domain of muon zenith angles
0095         idom = np.linspace(0,1, num_idom)   # domain of probability 0->1 which is domain of inverse cdf
0096         icdf = np.interp( idom, cdf, dom )  # interpolated invertion y<->x of the cdf
0097
0098         self.x = x        # domain of muon zenith angles
0099         self.y = y        # non-normalized PDF of muon zenith angles
0100         self.cdf = cdf
0101         self.icdf = icdf
0102         self.idom = idom
0103
0104
0105 class MuonAzimuthAngle(Generator):
0106     def __init__(self, x, num_idom=1000):
0107
0108         y = np.ones(len(x),dtype=np.float32)
0109         cdf = np.cumsum(y)
0110         cdf /= cdf[-1]       # CDFs tend to 1.
0111
0112         dom = x                             # domain of muon zenith angles
0113         idom = np.linspace(0,1, num_idom)   # domain of probability 0->1 which is domain of inverse cdf
0114         icdf = np.interp( idom, cdf, dom )  # interpolated invertion y<->x of the cdf
0115
0116         self.x = x
0117         self.y = y
0118         self.cdf = cdf
0119         self.icdf = icdf
0120         self.idom = idom
0121
0122
0123
0124 def plot3d(pos):
0125     pl = pv.Plotter()
0126     pl.add_points(pos)
0127     pl.show_grid()
0128     cp = pl.show()
0129     return cp
0130
0131
0132 if __name__ == '__main__':
0133
0134     zdom = np.linspace(0, 0.5*np.pi, 1000)
0135     adom = np.linspace(0, 2.*np.pi,  1000)
0136
0137     zen = MuonZenithAngle(zdom)
0138     azi = MuonAzimuthAngle(adom)
0139
0140     size = 10000
0141     u0 = np.random.uniform(size=size)
0142     u1 = np.random.uniform(size=size)
0143
0144     zenith = zen(u0)   # generate a sample of muon zenith angles that follows the desired distribution
0145     azimuth = azi(u1)
0146
0147     pos = np.zeros( (len(zenith),3), dtype=np.float32 )
0148
0149     pos[:,0] = np.sin(zenith)*np.cos(azimuth)
0150     pos[:,1] = np.sin(zenith)*np.sin(azimuth)
0151     pos[:,2] = np.cos(zenith)
0152
0153     plot3d(pos)
0154
0155     fig, axs = plt.subplots(2, 2)
0156
0157     zen.plot(axs[0][0])
0158     azi.plot(axs[0][1])
0159
0160     axs[1][0].hist(zenith, bins=500)
0161     axs[1][1].hist(azimuth, bins=500)
0162
0163     plt.show()
0164
0165     #pos = uniform_points_on_sphere(size=10_000)
0166
0167
0168
0169
0170
0171
0172
0173