Emission profiles


In [ ]:


In [1]:
import pschitt.emission as em
import pschitt.vizualisation as viz
import pschitt.sky_objects as sky
import pschitt.geometry as geo
import pschitt.camera_image as ci
import numpy as np
import matplotlib.pyplot as plt
import math

In [2]:
%matplotlib inline

In [ ]:

Here are the angular emission profiles accessible in the emission module.
These profiles may be used as a particle emission angular distribution.


In [13]:
fig, ax1 = plt.subplots(figsize=(10,6))
ax1 = viz.plot_angular_emission_profile(em.angular_profile.constant, 
                                        angle_max=2, label='constant', ax=ax1, color='black')
ax1 = viz.plot_angular_emission_profile(em.angular_profile.heaviside, 1,
                                        label='heaviside', ax=ax1)
ax1 = viz.plot_angular_emission_profile(em.angular_profile.exp_peak, 1.4, 10,
                                        label='exp_peak', ax=ax1)
ax1 = viz.plot_angular_emission_profile(em.angular_profile.lgdt06, 0.6,
                                        label='lgdt06', ax=ax1)
ax1.legend()


Out[13]:
<matplotlib.legend.Legend at 0x1a21def978>

In [ ]:

Demonstration on a real exemple

Setting up


In [15]:
shower = sky.shower()

shower.impact_point = np.array([50,0,0])

#shower direction
shower.alt = math.radians(70)
shower.az = math.radians(0)

#pointing direction off-axis
talt = math.radians(70)
taz = math.radians(0)

#shower parameters
shower_top = 12000
shower_length = 8000
shower_width = 100

shower.number_of_particles = int(1e4)

In [16]:
shower.gaussian_ellipsoide(shower_top, shower_length, shower_width)

In [17]:
tel_normal = geo.altaz_to_normal(talt,taz)
tel1 = geo.Telescope([0,0,0],tel_normal)
tel2 = geo.Telescope([100,0,0],tel_normal)
tel3 = geo.Telescope([200,0,0],tel_normal)
tel4 = geo.Telescope([400,0,0],tel_normal)
tel5 = geo.Telescope([0,100,0],tel_normal)

alltel = [tel1, tel2, tel3, tel4, tel5]

In [18]:
# alltel = geo.load_telescopes('../share/arrays/CTA_paranal.txt', normal=geo.altaz_to_normal(talt, taz))

Site Map


In [19]:
ax = viz.plot_array(alltel)
ax.scatter(shower.impact_point[0], shower.impact_point[1], color='black', label='Impact point', marker="+")
ax.legend()


Out[19]:
<matplotlib.legend.Legend at 0x1a21f7bf28>

In [ ]:

Producing images in telescopes


In [20]:
trigger_intensity = 20.
noise = 0

Constant profile


In [21]:
shower.particles_angular_emission_profile = em.angular_profile.constant
shower.particles_angular_emission_profile_kwargs = {}

In [22]:
ci.array_shower_imaging(shower, alltel, noise)

In [23]:
fig, axes = plt.subplots(1, len(alltel), figsize=(20,3))
for tel, ax in zip(alltel, axes):
    ax = viz.display_camera_image(tel, s=4, ax=ax);
    ax.set_title("Signal sum = {}".format(tel.signal_hist.sum()))



In [ ]:

Heaviside profile


In [24]:
shower.particles_angular_emission_profile = em.angular_profile.heaviside
shower.particles_angular_emission_profile_kwargs = {'limit':0.018}

In [25]:
ci.array_shower_imaging(shower, alltel, noise)

In [26]:
fig, axes = plt.subplots(1, len(alltel), figsize=(20,3))
for tel, ax in zip(alltel, axes):
    ax = viz.display_camera_image(tel, s=4, ax=ax);
    ax.set_title("Signal sum = {}".format(tel.signal_hist.sum()))


Only the telescope 5, right under the shower, got signal


In [ ]:

Exp peak profile


In [42]:
shower.particles_angular_emission_profile = em.angular_profile.exp_peak
shower.particles_angular_emission_profile_kwargs = {'eta':0.02, 'alpha':10}

In [43]:
ci.array_shower_imaging(shower, alltel, noise)

In [44]:
fig, axes = plt.subplots(1, len(alltel), figsize=(20,3))
for tel, ax in zip(alltel, axes):
    ax = viz.display_camera_image(tel, s=4, ax=ax);
    ax.set_title("Signal sum = {}".format(tel.signal_hist.sum()))


Maximum emission for telescopes at ~200m of the impact point


In [ ]:


In [ ]:

lgdt [Lemoine et al 06] profile


In [30]:
shower.particles_angular_emission_profile = em.angular_profile.lgdt06
shower.particles_angular_emission_profile_kwargs = {'eta':0.001}

In [31]:
ci.array_shower_imaging(shower, alltel, noise)

In [32]:
fig, axes = plt.subplots(1, len(alltel), figsize=(20,3))
for tel, ax in zip(alltel, axes):
    ax = viz.display_camera_image(tel, s=4, ax=ax);
    ax.set_title("Signal sum = {}".format(tel.signal_hist.sum()))


Telescopes outside the cone receive less signal


In [ ]:


In [ ]: