Number of event

import numpy as np
import matplotlib.pyplot as plt

def asymmetry():
    B_values = np.linspace(0.00, 0.4,51)  # an array of B values from 0.01 to 1.00
    N_event_array = []
    A_event_array = []
    for B in B_values:
        N_point_array = []
        A_point_array = []
        Nplus_array = np.arange(1000, 10000, 1) # an array of Nplus values from 10 to 1000
        Nminus_array = (1 - B) * Nplus_array
        N_array = np.round(Nplus_array + Nminus_array,2)  # total number
        A_array = np.round((Nplus_array - Nminus_array) / (Nplus_array + Nminus_array),5) # calculate the value of asymmetry A
        
        # Calculte the uncertainty of Nplus and Nminus
        sigma_Nplus_array = np.sqrt(Nplus_array) 
        sigma_Nminus_array = np.sqrt(1 - B) * sigma_Nplus_array

        # Calculate the uncertainty of A due to the uncertainty of N, where dA/dNplus and dA/dNminus are the partial derivative
        #of A with respect to Nplus and Nminus and then add these contribution to the uncertainty of A

        sigma_A1_array = ((2 * Nminus_array * sigma_Nplus_array) / (Nplus_array + Nminus_array)**2)**2 #using the formula (dA/dNplus)^2*(sigmaNplus)^2
        sigma_A2_array = (-(2 * Nplus_array * sigma_Nminus_array) / (Nplus_array + Nminus_array)**2)**2 #using the formula (dA/dNminus)^2*(sigmaNplus)^2

        delta_A_array = np.round(np.sqrt(sigma_A1_array + sigma_A2_array),5) # the uncertainty of A


        min_index = np.argmin(delta_A_array)
       
        max_index = np.argmax(delta_A_array)
       
        n = 2e-5
        if 0.0<B<=0.06:
            n=n
        elif 0.06<B<=0.14:
            n=n*10
        elif 0.14<B<=0.36:
            n=n*1e2
        elif 0.36<B<=0.7:
            n=n*3e3
        elif 0.7<B<0.82:
            n=n*1e4
        
        #Find the number of event at intersection between Asymmetry and Uncertainty for every value of B, and save them to new array
        for i in range(len(A_array)):
            if A_array[i] < np.round(delta_A_array[min_index],5) or A_array[i] > delta_A_array[max_index]:
                N_point_array.append(0)
            elif abs(A_array[i] - delta_A_array[i]) <=n:
                N_point = N_array[i]
                A_point = A_array[i]
                N_point_array.append(N_point)
                A_point_array.append(A_point)
        
        if len(N_point_array) == 0:
            N_point_array.append(0)
        if len(A_point_array) == 0:
            A_point_array.append(0)
        
        A_event = A_point_array[0]
        N_event = N_point_array[0]
        for i in range(1, len(N_point_array)):
            if N_point_array[i] > N_event:
                N_event = N_point_array[i]
                A_event = A_point_array[i]
        
        N_event_array.append(N_event)
        A_event_array.append(A_event)
       
    
    while 0 in N_event_array:
        N_event_array.remove(0)
        A_event_array.remove(0)
    
    fig, ax = plt.subplots()
    ax.plot(N_event_array, A_event_array, label='Uncertainty')
    ax.set_xlabel('Number of event')
    ax.set_ylabel('Asymmetry and Uncertainty')
    ax.set_title(f'Asymmetry and Uncertainty vs. N_array for B')
    ax.legend()
    plt.show()
asymmetry()