Python Params.matFunction方法代码示例

# 需要导入模块: import Params [as 别名]
# 或者: from Params import matFunction [as 别名]
def getScatterAnglePDG(x, dt):
    # return thetax, thetay, yx, yy
    # given a velocity and timestep, compute a
    # deflection angle and deviation due to multiple scattering
    # Update the position/velocity and return deflection
    #
    # x is a 6-element vector (x,y,z,px,py,pz)
    # returns deflection angles thetax, thetay and
    # displacements yx, yy in two orthogonal directions to p

    p = x[3:]
    magp = np.linalg.norm(p) # must be in MeV
    E = np.sqrt(magp**2 + Params.m**2)
    v = p/E
    beta = np.linalg.norm(v)

    dx = (beta*2.9979e-1) * dt  # in m
    
    ## in the following we generate a random deflection angle theta
    ## and transverse displacement y for the given momentum. This is taken
    ## from the PDG review chapter on the Passage of Particles through Matter

    mat = Params.matFunction(x[0],x[1],x[2])

    X0 = Params.materials[mat][3]

    if X0<=0:
        return np.zeros(6)

    # rms of projected theta distribution.
    theta0 = 13.6/(beta*magp) * abs(Params.Q) * np.sqrt(dx/X0) * (1 + 0.038*np.log(dx/X0))
    
    # correlation coefficient between theta_plane and y_plane
    rho = 0.87
    
    getRandom = np.random.normal

    z1 = getRandom()
    z2 = getRandom()
    yx = z1*dx*theta0 * np.sqrt((1-rho**2)/3) + z2*rho*dx*theta0/np.sqrt(3)
    thetax = z2*theta0
    
    z1 = getRandom()
    z2 = getRandom()
    yy = z1*dx*theta0 * np.sqrt((1-rho**2)/3) + z2*rho*dx*theta0/np.sqrt(3)
    thetay = z2*theta0

    return thetax, thetay, yx, yy

# 需要导入模块: import Params [as 别名]
# 或者: from Params import matFunction [as 别名]
def multipleScatterKuhn(x, dt):
    # use the method from Kuhn paper

    if Params.matFunction(x[0],x[1],x[2])=='air':
        return np.zeros(6)

    p = x[3:]
    theta = getScatterAngleKuhn(x, dt)

    if theta==-1:
        return multipleScatterPDG(x,dt)

    vx = getNormVector(p)
    
    # deflection in momentum
    defl = np.linalg.norm(p) * (theta*vx)

    return np.append(np.zeros(3), defl)

# 需要导入模块: import Params [as 别名]
# 或者: from Params import matFunction [as 别名]
def multipleScatterPDG(x, dt):
    # get the angles/displacements from above function and return the
    # net change in x=(x,y,z,px,py,pz)

    if Params.matFunction(x[0],x[1],x[2])=='air':
        return np.zeros(6)

    p = x[3:]

    vx = getNormVector(p)
    vy = np.cross(vx, p/np.linalg.norm(p))

    thetax, thetay, yx, yy = getScatterAnglePDG(x, dt)

    # transverse displacement
    disp = yx*vx + yy*vy
    
    # deflection in momentum
    defl = np.linalg.norm(p) * (thetax*vx + thetay*vy)

    return np.append(disp, defl)

# 需要导入模块: import Params [as 别名]
# 或者: from Params import matFunction [as 别名]
def getKuhnScatteringParams(x, dt):

    mat = Params.matFunction(x[0],x[1],x[2])

    Z,A,rho,X0 = Params.materials[mat]

    z = abs(Params.Q)

    p = x[3:]
    magp = np.linalg.norm(p)
    v = p/np.sqrt(magp**2 + Params.m**2)
    beta = np.linalg.norm(v)

    ds = beta * 2.9979e1 * dt

    Xc = np.sqrt(0.1569 * z**2 * Z*(Z+1) * rho * ds / (magp**2 * beta**2 * A))
    b = np.log(6700*z**2*Z**(1./3)*(Z+1)*rho*ds/A / (beta**2+1.77e-4*z**2*Z**2))

    if b<3:
        if not Params.MSCWarning:
            print "Warning: something (probably Q) is too small! Using PDG MSC algorithm."
            Params.MSCWarning = True
        return -1,-1

    ## we want to solve the equation B-log(B) = b. Using Newton-Raphson

    B = b
    prevB = 2*B
    
    f = lambda x: x-np.log(x)-b
    fp = lambda x: 1-1./x

    while abs((B-prevB)/prevB)>0.001:
        prevB = B
        B = B - f(B)/fp(B)

        
    # use B+1 for correction at intermediate angles
    return Xc, B+1

# 需要导入模块: import Params [as 别名]
# 或者: from Params import matFunction [as 别名]
def doEnergyLoss(x, dt):
    ## returns new x after losing proper amount of energy according to Bethe-Bloch

    p = x[3:]
    magp = np.linalg.norm(p)
    E = np.sqrt(magp**2+Params.m**2)
    gamma = E/Params.m
    beta = magp/E;
    me = 0.511;  #electron mass in MeV
    
    Wmax = 2*me*beta**2*gamma**2/(1+2*gamma*me/Params.m + (me/Params.m)**2)
    K = 0.307075  # in MeV cm^2/mol

    mat = Params.matFunction(x[0],x[1],x[2])
    Z,A,rho,X0 = Params.materials[mat]
    I,a,k,x0,x1,Cbar,delta0 = Params.dEdx_params[mat]

    I = I/1e6  ## convert from eV to MeV

    xp = np.log10(magp/Params.m)
    if xp>=x1:
        delta = 2*np.log(10)*xp - Cbar
    elif xp>=x0:
        delta = 2*np.log(10)*xp - Cbar + a*(x1-xp)**k
    else:
        delta = delta0*10**(2*(xp-x0))

    # mean energy loss in MeV/cm
    dEdx = K*rho*Params.Q**2*Z/A/beta**2*(0.5*np.log(2*me*beta**2*gamma**2*Wmax/I**2) - beta**2 - delta/2)

    dE = dEdx * beta*2.9979e1 * dt

    if dE>(E-Params.m):
        return np.array([x[0], x[1], x[2], 0, 0, 0])

    newmagp = np.sqrt((E-dE)**2-Params.m**2)
    x[3:] = p*newmagp/magp

    return x