<jupyter_start><jupyter_text>Aeff toy MC<jupyter_code>import numpy as np
from matplotlib import pyplot as plt

plt.rcParams['text.usetex'] = False<jupyter_output><empty_output><jupyter_text>define the function for the selection efficiency, as a function of log10(E):<jupyter_code>def Eeff(logE):
    x0=2.7
    y=1-1./(np.exp((logE-x0)/0.5)+1)  
    
    return y<jupyter_output><empty_output><jupyter_text>generatig random distributions using numpy arrays:<jupyter_code>Nev=10000000 # n. of events we simulate
E_MIN=100
E_MAX=10000000

# generate energies with E^-1 spectrum, using the inversion method 
u1=np.random.rand(Nev)
energies=E_MIN*np.exp(u1*np.log(E_MAX/E_MIN))

#generate phi angles (even if this is not used at the moment)
u2=np.random.rand(Nev)
phi = 2*np.pi*u2

#generate theta angles 
u3=np.random.rand(Nev)
theta= np.arccos(u3)

#generate x,y impact points with -1<x<1, -1<y<1  (2x2 m^2) 
x= -1+2*np.random.rand(Nev)
y=-1+2*np.random.rand(Nev)
     
 
#theta1=[40/180 *np.pi for i in range(0,int(0.6*Nev))]
#theta2=[60/180 *np.pi for i in range(0,int(0.4*Nev))]
#theta=np.concatenate((theta1,theta2))<jupyter_output><empty_output><jupyter_text>apply cuts using array masks:<jupyter_code>#energy mask using the hit/miss method. 
u5=np.random.rand(Nev)
pe=Eeff(np.log10(energies))
e_mask=u5<pe # keep the cases where u < y 

#combine all cuts:
totalMask=(e_mask& (x>-np.sqrt(2)/2)& (x<np.sqrt(2)/2)  & (y>-np.sqrt(2)/2)& (y<np.sqrt(2)/2))# detector dimensions -0.5-0.5  on both axis<jupyter_output><empty_output><jupyter_text>plot the distributions:Energy: (note that the resulting histo include all cuts, not only on energy)<jupyter_code>plt.figure(1,figsize=(7,5))
plt.hist(np.log10(energies[totalMask]), bins=50, alpha=1, histtype='step', label='final E')
plt.hist(np.log10(energies), bins=50, alpha=1,  histtype='step',label='simulated E')
plt.legend()
plt.xlabel('log10(E)')
plt.ylabel('counts')
plt.savefig("counts(E).png")
plt.show()<jupyter_output><empty_output><jupyter_text>theta and phi:<jupyter_code>plt.figure(2,figsize=(7,5))
plt.hist(theta[totalMask], bins=50, alpha=1, histtype='step', label='theta final')
plt.hist(theta, bins=50, alpha=1, histtype='step', label='theta simulated')
#plt.figure(3)
plt.hist(phi[totalMask], bins=50, alpha=1, histtype='step', label='phi final')
plt.hist(phi, bins=50, alpha=1, histtype='step', label='phi simulated')
plt.xlabel('[rad]')
plt.ylabel('counts')
plt.legend()
plt.savefig("counts(angle).png")
plt.show()<jupyter_output><empty_output><jupyter_text>x,y<jupyter_code>plt.figure(3,figsize=(7,5))
plt.hist(x[totalMask], bins=50, alpha=1, histtype='step', label='x final')
plt.hist(x, bins=50, alpha=1, histtype='step', label='x simulated')
#plt.figure(3)
plt.hist(y[totalMask], bins=50, alpha=1, histtype='step', label='y final')
plt.hist(y, bins=50, alpha=1, histtype='step', label='y simulated')
plt.xlabel('[m]')
plt.ylabel('counts')
plt.legend()
plt.savefig("counts(area).png")
plt.show()<jupyter_output><empty_output><jupyter_text>create rensponse matrix:<jupyter_code>nBins2d=40
H,xbins,ybins=np.histogram2d(np.log10(energies[totalMask]),np.cos(theta[totalMask]),nBins2d)

# histogram normaliztion
normFactor=Nev/(2*nBins2d**2) # the factor 4 comes form the area reatio
H=H/(normFactor)
# plot
#y=[i/np.pi*180 for i in ybins]
fig, ax1 = plt.subplots(figsize=(8, 5))
p=ax1.pcolormesh(xbins,ybins,H.T, cmap='rainbow')
bar=fig.colorbar(p, ax=ax1, extend='both')
bar.ax.set_ylabel('Aeff [m^2]', rotation=270,labelpad=13,fontsize=13)
ax1.set_xlabel('log10(E)',fontsize=14)
ax1.set_ylabel('theta[°]',fontsize=14)
plt.savefig("Aeff_costheta.png")
plt.show()
y
ybins
nBins2d=40
H,xbins,ybins=np.histogram2d(np.log10(energies[totalMask]),theta[totalMask],nBins2d)

# histogram normaliztion
normFactor=Nev/(2*nBins2d**2) # the factor 4 comes form the area reatio
H=H/(normFactor)
# plot
y=[i/np.pi*180 for i in ybins]
fig, ax1 = plt.subplots(figsize=(8, 5))
p=ax1.pcolormesh(xbins,y,H.T, cmap='rainbow')
bar=fig.colorbar(p, ax=ax1, extend='both')
bar.ax.set_ylabel('Aeff [m^2]', rotation=270,labelpad=13,fontsize=13)
ax1.set_xlabel('log10(E)',fontsize=14)
ax1.set_ylabel('theta[°]',fontsize=14)
#ax1.set_xscale("log")
plt.savefig("Aeff_theta.png")
plt.show()

# Finallly slect one bin in cos(theta) and plot Aeff vs log10(E) for that bin

plt.figure(5)
halfxbin=(xbins[1]-xbins[0])/2.
plt.plot(xbins[:-1]+halfxbin, H.T[18],label="40°" )
plt.plot(xbins[:-1]+halfxbin, H.T[27],label="60°" )
plt.legend()
plt.xlabel('log10(E)')
plt.ylabel('Aeff [m]')
plt.savefig("40_60")
plt.show()

mean=[0.6*H.T[18][i]+0.4*H.T[27][i] for i in range(0,len(H.T[18]))]
# Finallly slect one bin in cos(theta) and plot Aeff vs log10(E) for that bin

plt.figure(5)
halfxbin=(xbins[1]-xbins[0])/2.
plt.plot(xbins[:-1]+halfxbin, mean )
plt.title("Mean" )
plt.xlabel('log10(E)')
plt.ylabel('Aeff [m]')
plt.savefig("Mean")
plt.show()

# Finallly slect one bin in cos(theta) and plot Aeff vs log10(E) for that bin

plt.figure(5)
halfxbin=(xbins[1]-xbins[0])/2.
plt.plot(xbins[:-1]+halfxbin, H.T[18],label="40°" )
plt.plot(xbins[:-1]+halfxbin, H.T[27],label="60°" )
plt.plot(xbins[:-1]+halfxbin,mean,label="mean" )
plt.legend()
plt.xlabel('log10(E)')
plt.ylabel('Aeff [m]')
plt.savefig("Aeff_all")
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.195260663507109, 'fraction_numerical': 0.05004739336492891, 'mean_word_length': 7.348101265822785, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '10245', 'license': 'MIT', 'n_tokens_mistral': 2482, 'n_tokens_neox': 2297, 'n_words': 381, 'path': 'Data_Analysis/Fenu_Maldera/Esercizio_1.ipynb', 'repo_name': 'andreasemeraro/MPM_Space_Sciences', 'size': '5275'}