Commit 6ed94a21 authored by Alice Donini's avatar Alice Donini
Browse files

update notebooks

parent 8419c545
Loading
Loading
Loading
Loading
+564 −0
Original line number Diff line number Diff line
%% Cell type:markdown id:ultimate-brazilian tags:

# Fast analysis from DL2

%% Cell type:markdown id:standing-cancellation tags:

This notebook explores DL2 **source independent** data and allows to do a fast analysis of the data.

%% Cell type:markdown id:outdoor-ridge tags:

## Import necessary packages:

%% Cell type:code id:smaller-marker tags:

``` python
%matplotlib inline
import os
import pandas as pd
import numpy as np
import re
import astropy.units as u
import matplotlib.pyplot as plt
import matplotlib.style as style
from matplotlib.offsetbox import AnchoredText

from lstchain.io.io import dl2_params_lstcam_key
from lstchain.reco.utils import get_effective_time, radec_to_camera
from lstchain.reco.utils import compute_theta2, extract_source_position, clip_alt
from ctapipe.containers import EventType
from ctapipe.coordinates import CameraFrame
from astropy.coordinates import AltAz, SkyCoord, EarthLocation
from astropy.coordinates import ICRS, Galactic, FK4, FK5
from astropy.time import Time
from gammapy.stats import WStatCountsStatistic

style.use('tableau-colorblind10')
plt.rcParams['font.size'] = 20
```

%% Cell type:code id:ba0c0523-4c9e-4031-bfc5-e61fb8d6f671 tags:

``` python
import lstchain
import gammapy
print("lstchain:", lstchain.__version__)
print("gammapy:", gammapy.__version__)
```

%% Cell type:markdown id:heavy-lunch tags:

## This function gets the source position in camera coordinates for any source

%% Cell type:markdown id:unknown-humanitarian tags:

Alternatively you can use lstchain.reco.extract_source_position() just giving the source name (p.e. "Crab")

%% Cell type:code id:saved-suite tags:

``` python
def extract_source_position_from_coord(
    data, coord, equivalent_focal_length=28 * u.m
):
    """
    Extract source position from data
    Parameters:
    -----------
    pandas.DataFrame data: input data
    str observed_source_name: Name of the observed source
    astropy.units.m equivalent_focal_length: Equivalent focal length of a telescope
    Returns:
    --------
    2D array of coordinates of the source in form [(x),(y)] in astropy.units.m
    """

    obstime = pd.to_datetime(data["dragon_time"], unit="s")
    pointing_alt = u.Quantity(data["alt_tel"], u.rad, copy=False)
    pointing_az = u.Quantity(data["az_tel"], u.rad, copy=False)
    source_pos_camera = radec_to_camera(
        coord,
        obstime,
        pointing_alt,
        pointing_az,
        focal=equivalent_focal_length,
    )
    source_position = [source_pos_camera.x, source_pos_camera.y]
    return source_position
```

%% Cell type:markdown id:foster-decision tags:

## Choose your input files

%% Cell type:markdown id:bc6fcb96 tags:

### If you use single data DL2 file:

%% Cell type:code id:7976824c tags:

``` python
#input_file='../data/src_indep/dl2_LST-1.Run2968_first10.h5'
#data=pd.read_hdf(input_file, key=dl2_params_lstcam_key)
```

%% Cell type:markdown id:e5adc197 tags:

### If you want to test this notebook with several full runs

%% Cell type:code id:72baa6c0 tags:

``` python
base_dir = '/fefs/aswg/workspace/alice.donini/Analysis/data/DL2/GRB221016A/20221016/v0.9/tailcut84/'
!ls $base_dir
```

%% Cell type:code id:dangerous-continuity tags:

``` python
%time
data=pd.DataFrame()
#runs=['8501']
#runs=['8863','8864','8865','8866','8867','8878', '8879','8880','8881'] # You can concatenate several runs

# find run number automatically
runs = []
for fname in os.listdir(base_dir):
    res = re.findall("dl2_LST-1.Run0(\d+).h5", fname)
    #print(res[0])
    if not res: continue
    runs.append(res[0])
runs.sort()
print(runs)
```

%% Cell type:code id:eb0245f4-57b4-4b62-a9f3-8ddc6455444c tags:

``` python
for run in runs:
    input_file = base_dir + 'dl2_LST-1.Run0' + run + '.h5'
    data=pd.concat([data, pd.read_hdf(input_file, key=dl2_params_lstcam_key)])
```

%% Cell type:markdown id:3a83fe7b tags:

## Take a look at the dataframe

%% Cell type:code id:01b3ae1c tags:

``` python
data
```

%% Cell type:code id:d8c6998d tags:

``` python
data.keys()
```

%% Cell type:markdown id:infinite-demonstration tags:

## Get the observation time and store some useful variables

%% Cell type:code id:arabic-testing tags:

``` python
obstime_real = get_effective_time(data)[0]
gammaness = np.array(data.gammaness)
leakage_intensity_width_2 = np.array(data.leakage_intensity_width_2)
intensity = np.array(data.intensity)
alt = np.array(data.alt_tel)
wl = np.array(data.wl)
event_type = np.array(data.event_type)
event_id = np.array(data.event_id)
obstime_real.to(u.min), obstime_real.to(u.h)
```

%% Cell type:markdown id:sized-lesson tags:

## Define some cuts

%% Cell type:code id:pressing-liechtenstein tags:

``` python
gammaness_cut=0.6
intensity_cut=50

intensity_cut_high=1e10
alt_min=0 * np.pi / 180
wl_cut=0.0
THETA2_GLOBAL_CUT=0.04
theta2_range=(0,1)
norm_range_theta2_min=0.2
norm_range_theta2_max=0.5
```

%% Cell type:code id:opposite-brother tags:

``` python
condition = (gammaness > gammaness_cut) \
                                     & (intensity > intensity_cut) \
                                     & (intensity < intensity_cut_high) \
                                     & (alt > alt_min) \
                                     & (wl > wl_cut) \
                                     & (event_type != EventType.FLATFIELD.value)\
                                     & (event_type != EventType.SKY_PEDESTAL.value)\
                                     & (leakage_intensity_width_2 < 0.2)
```

%% Cell type:code id:aboriginal-convenience tags:

``` python
selected_data=data[condition] #We apply now the cuts so the calculations will be quicker
```

%% Cell type:markdown id:happy-monaco tags:

## Define the source coordinates

%% Cell type:code id:polyphonic-complexity tags:

``` python
#coords = ["05 34 31.94 +22 00 52.2"] #Crab coordinates
#coords = ["14 05 39.268  +20 38 28.29"]
#coordinates = SkyCoord(coords, frame=ICRS, unit=(u.hourangle, u.deg))

# Ra, Dec in deg
source_name = 'GRB221016A'
coords = [38.9491, -34.6241]
coordinates = SkyCoord(ra=coords[0],dec=coords[1], frame=ICRS, unit=u.deg)

filename_output = f'{source_name}'
coordinates
```

%% Cell type:markdown id:planned-fence tags:

## Calculate the source true position

%% Cell type:code id:compressed-vertical tags:

``` python
%%time
true_source_position = extract_source_position_from_coord(selected_data, coordinates)
#true_source_position = extract_source_position(selected_data, 'Crab') #If its a catalogued source, like the Crab, you can use this lstchain function
off_source_position = [element * -1 for element in true_source_position]
```

%% Cell type:markdown id:usual-physics tags:

## Compute theta2 of the ON and OFF data

%% Cell type:code id:shaped-exhibit tags:

``` python
theta2_on = np.array(compute_theta2(selected_data, true_source_position))
theta2_off = np.array(compute_theta2(selected_data, off_source_position))
```

%% Cell type:markdown id:secure-supplier tags:

## Create histograms for theta2 plots

%% Cell type:code id:vanilla-killer tags:

``` python
#nbins=100 #Choose your preferred number of bins
nbins=round((theta2_range[1]/THETA2_GLOBAL_CUT)*2) # Make the histogram so there are only two bins before the theta2 cut
hist_on, bin_edges_on=np.histogram(theta2_on,density=False, bins=nbins, range=theta2_range)
hist_off, bin_edges_off=np.histogram(theta2_off, density=False, bins=nbins, range=theta2_range)

bin_width=bin_edges_on[1]-bin_edges_off[0]
bin_center=bin_edges_on[:-1]+(bin_width/2)
```

%% Cell type:code id:31a5c1fa tags:

``` python
#plt.hist(hist_on)
```

%% Cell type:markdown id:forward-anthony tags:

## Calculate the Li&Ma significance

%% Cell type:code id:b8983e58 tags:

``` python
N_on = np.sum(hist_on[bin_edges_on[1:]<=THETA2_GLOBAL_CUT])
N_off = np.sum(hist_off[bin_edges_off[1:]<=THETA2_GLOBAL_CUT])
N_off_err = np.sqrt(N_off)

N_on,N_off,N_off_err
```

%% Cell type:code id:e2ee1e28 tags:

``` python
# argmin restituisce l'indice del minimo valore dell'array
idx_min = (np.abs(bin_edges_on - norm_range_theta2_min)).argmin()
idx_max = (np.abs(bin_edges_on - norm_range_theta2_max)).argmin()
idx_min,idx_max
```

%% Cell type:code id:b1cba83a tags:

``` python
Non_norm = np.sum(hist_on[idx_min:idx_max])
Noff_norm = np.sum(hist_off[idx_min:idx_max])
alpha = Non_norm/Noff_norm
alpha, 1/alpha
```

%% Cell type:code id:f5dad613 tags:

``` python
#plt.hist(hist_off[idx_min:idx_max])
```

%% Cell type:code id:lyric-desperate tags:

``` python
stat = WStatCountsStatistic(n_on=N_on, n_off=N_off, alpha=alpha)
significance_lima = stat.sqrt_ts

textstr = r'N$_{{\rm on}}$ = {:.0f} '\
            f'\n'\
            r'N$_{{\rm off}}$ = {:.0f} $\pm$ {:.0f}'\
            f'\n'\
            r'Time = {:.1f}'\
            f'\n'\
            r'1/$\alpha$ = {:.2f}'\
            f'\n'\
            r'LiMa Significance = {:.1f} $\sigma$ '.format(N_on,
                                                      N_off,
                                                      N_off_err,
                                                      obstime_real.to(u.h),
                                                      1/alpha,
                                                      significance_lima)
```

%% Cell type:markdown id:tired-ottawa tags:

## Plot Theta2

%% Cell type:code id:stopped-sigma tags:

``` python
fig, ax = plt.subplots(figsize=(12, 10))

ax.errorbar(bin_center, hist_on, yerr=np.sqrt(hist_on), fmt='o', label='ON data', ms=10, color='tab:orange')
ax.errorbar(bin_center, hist_off, yerr=np.sqrt(hist_off),fmt='s',label='Background', ms=10, color='tab:blue')
ax.set_xlim(0, 0.5)
#ax.set_xlim(left=0)
ax.grid(ls='dashed')
ax.axvline(THETA2_GLOBAL_CUT, color='black',ls='--',alpha=0.75)
ax.set_xlabel("$\\theta^{2} [deg^{2}]$")
ax.set_ylabel("Counts")
ax.legend(loc='lower right')
ax.set_title(source_name)

box_prop = dict(boxstyle='Round', facecolor='white', alpha=0.5)
text_prop = dict(fontsize=18, bbox=box_prop)
txt = AnchoredText(textstr, loc=1, transform=ax.transAxes, prop=text_prop, frameon=False)
ax.add_artist(txt)

plt.savefig(os.path.join(base_dir, f'{filename_output}_gammaness{gammaness_cut}_intensity{intensity_cut}_fast_theta2.png'), dpi=300)
```

%% Cell type:markdown id:skilled-attitude tags:

## Plot excess counts

%% Cell type:code id:second-consideration tags:

``` python
excess=hist_on-hist_off
excess_err=np.sqrt(np.sqrt(excess**2))
```

%% Cell type:code id:mechanical-norwegian tags:

``` python
fig, ax = plt.subplots(figsize=(12, 10))

ax.errorbar(bin_center, excess, yerr=excess_err,fmt='o',color='tab:blue',label='Excess counts')
ax.bar(bin_edges_on[:-1], excess, width = bin_width, align='edge', color='lightblue', ec='lightblue', alpha=0.5)
ax.axhline(0, color='darkgray')
ax.set_xlim(0, 0.5)
#ax.set_xlim(left=0)
ax.grid()
ax.axvline(THETA2_GLOBAL_CUT, color='black', ls='--', alpha=0.75)
ax.grid(ls='dashed')
ax.set_xlabel("$\\theta^{2} [deg^{2}]$")
ax.set_ylabel("Counts")
ax.legend(title=f'Significance = {significance_lima:.1f} $\sigma$')

plt.savefig(os.path.join(base_dir, f'{filename_output}_gammaness{gammaness_cut}_intensity{intensity_cut}_fast_excess.png'), dpi=300)
```

%% Cell type:markdown id:86bfa438 tags:

# We can make theta2 plots for different bins in energy

%% Cell type:markdown id:718e7e80 tags:

## Define the energy binning

%% Cell type:code id:2cda29d5 tags:

``` python
log_reco_energy = np.array(selected_data.log_reco_energy)
emin=0.01 * u.TeV
emax=10 * u.TeV
n_bins_energy=3
log_energy = np.linspace(np.log10(emin.to_value()),
                         np.log10(emax.to_value()),
                         n_bins_energy + 1)
```

%% Cell type:markdown id:f1935cbf tags:

## Now calculate the significante and produce a theta2 plot for each energy bin

%% Cell type:code id:e8242e86 tags:

``` python
for i in range(n_bins_energy):
    condition_energy_bin= (log_reco_energy < log_energy[i+1]) \
        & (log_reco_energy >= log_energy[i])
    data_bin=selected_data[condition_energy_bin]
    theta2_on = np.array(compute_theta2(data_bin, (true_source_position[0][condition_energy_bin], true_source_position[1][condition_energy_bin])))
    theta2_off = np.array(compute_theta2(data_bin, (off_source_position[0][condition_energy_bin], off_source_position[1][condition_energy_bin])))

    hist_on, bin_edges_on=np.histogram(theta2_on,density=False, bins=nbins, range=theta2_range)
    hist_off, bin_edges_off=np.histogram(theta2_off, density=False, bins=nbins, range=theta2_range)

    bin_width=bin_edges_on[1]-bin_edges_off[0]
    bin_center=bin_edges_on[:-1]+(bin_width/2)
    N_on = np.sum(hist_on[bin_edges_on[1:]<=THETA2_GLOBAL_CUT])
    N_off = np.sum(hist_off[bin_edges_off[1:]<=THETA2_GLOBAL_CUT])

    idx_min = (np.abs(bin_edges_on - norm_range_theta2_min)).argmin()
    idx_max = (np.abs(bin_edges_on - norm_range_theta2_max)).argmin()

    Non_norm = np.sum(hist_on[idx_min:idx_max])
    Noff_norm = np.sum(hist_off[idx_min:idx_max])

    alpha =  Non_norm/Noff_norm

    stat = WStatCountsStatistic(n_on=N_on, n_off=N_off, alpha=alpha)
    significance_lima = stat.sqrt_ts

    textstr = r'N$_{{\rm on}}$ = {:.0f} '\
            f'\n'\
            r'N$_{{\rm off}}$ = {:.0f} $\pm$ {:.0f}'\
            f'\n'\
            r'Time = {:.1f}'\
            f'\n'\
            r'LiMa Significance = {:.1f} $\sigma$ '.format(N_on,
                                                      N_off,
                                                      N_off_err,
                                                      obstime_real.to(u.h),
                                                      significance_lima)

    props = dict(boxstyle='round', facecolor='wheat', alpha=0.95)

    fig, ax = plt.subplots(figsize=(12, 10))

    ax.errorbar(bin_center, hist_on, yerr=np.sqrt(hist_on), fmt='o', label='ON data', ms=10, color='tab:orange')
    ax.errorbar(bin_center, hist_off, yerr=np.sqrt(hist_off),fmt='s',label='Background', ms=10, color='tab:blue')
    ax.set_xlim(0, 0.5)
    ax.grid(ls='dashed')
    ax.axvline(THETA2_GLOBAL_CUT, color='black',ls='--',alpha=0.75)
    ax.set_xlabel("$\\theta^{2} [deg^{2}]$")
    ax.set_ylabel("Counts")
    ax.legend(bbox_to_anchor=(0.99, 0.8))
    ax.set_title(r'{:.2f} TeV to {:.2f} TeV'.format(10**log_energy[i], 10**log_energy[i+1]))

    box_prop = dict(boxstyle='Round', facecolor='white', alpha=0.5)
    text_prop = dict(fontsize=20, bbox=box_prop)
    txt = AnchoredText(textstr, loc=1, transform=ax.transAxes, prop=text_prop, frameon=False)
    ax.add_artist(txt)
```

%% Cell type:markdown id:electoral-substitute tags:

## Plot the position of the selected events

%% Cell type:markdown id:featured-siemens tags:

### In camera coordinates

%% Cell type:code id:arbitrary-karma tags:

``` python
fig, ax = plt.subplots(figsize=(12, 12))

skymap=ax.hist2d(selected_data['reco_src_x'],
           selected_data['reco_src_y'], bins=100,
           range=[(-1,1), (-1,1)])
ax.set_xlabel("x (m)")
ax.set_ylabel("y (m)")
```

%% Cell type:markdown id:700caf5a tags:

You can see the wobble positions if you use one more than one run

%% Cell type:markdown id:friendly-mambo tags:

### In Sky coordinates

%% Cell type:code id:italian-basics tags:

``` python
location = EarthLocation.from_geodetic(-17.89139 * u.deg, 28.76139 * u.deg, 2184 * u.m) #Location of LST1
obstime = pd.to_datetime(selected_data["dragon_time"], unit="s")
horizon_frame = AltAz(location=location, obstime=obstime)
```

%% Cell type:code id:explicit-recorder tags:

``` python
%%time
pointing_alt = u.Quantity(selected_data["alt_tel"], u.rad, copy=False)
pointing_az = u.Quantity(selected_data["az_tel"], u.rad, copy=False)
pointing_direction=SkyCoord(alt=clip_alt(pointing_alt), az=pointing_az, frame=horizon_frame)
camera_frame = CameraFrame(focal_length=28 * u.m,
                           telescope_pointing=pointing_direction,
                           obstime=obstime,
                           location=location)
```

%% Cell type:code id:civil-baking tags:

``` python
camera_coords = SkyCoord(x=selected_data['reco_src_x'], y=selected_data['reco_src_y'], frame=camera_frame, unit=(u.m, u.m))
```

%% Cell type:code id:assumed-visibility tags:

``` python
radec_coords=camera_coords.transform_to(frame=ICRS)
```

%% Cell type:code id:9f8e330d tags:

``` python
coordinates.ra
```

%% Cell type:code id:horizontal-asbestos tags:

``` python
fig, ax = plt.subplots(figsize=(12, 12))
radec_skymap=plt.hist2d(radec_coords.ra.to_value(u.hourangle),
           radec_coords.dec.value,
           bins=100, cmap='hot')
ax.plot(coordinates.ra.to_value(u.hourangle), coordinates.dec.value, marker='+',
        color='lime', markersize=30)
ax.text(coordinates.ra.to_value(u.hourangle)+0.01, coordinates.dec.value+0.1, 'Source', color='lime')
ax.grid(color='white')
ax.set_xlabel("Right Ascension (h)")
ax.set_ylabel("Declination (deg)")
```

notebooks/Fast_analysis.ipynb

deleted100644 → 0
+0 −1273

File deleted.

Preview size limit exceeded, changes collapsed.

+1 −2

File changed.

Preview size limit exceeded, changes collapsed.