preprocess_hdf5.py

from ctapipe.calib.camera.r1 import HessioR1Calibrator
from ctapipe.calib.camera.dl0 import CameraDL0Reducer
from ctapipe.calib.camera import CameraDL1Calibrator
from ctapipe.image.charge_extractors import LocalPeakIntegrator
from ctapipe.image import tailcuts_clean, hillas_parameters
from ctapipe.io.hessio import hessio_event_source
import pickle
import os
from multiprocessing import Pool
import sys
import pandas as pd
import numpy as np
import time
import h5py
from tqdm import tqdm


def string_len(string, length):
    string = str(string)
    ersetzung = " "
    if length == 2:
        ersetzung = "0"
    if len(string) < length:
        string = ersetzung + string
    return string


def get_num_events(Filename, right_tel):
    try:
        source = hessio_event_source(Filename, allowed_tels=right_tel)
    except:
        os.exit(1)

    num_events = 0
    for event in source:
        num_events += 1
    return num_events


def set_right_tel(Filename):
    right_tel = []
    source = hessio_event_source(Filename)
    for event in source:
        for tel_id in event.r0.tels_with_data:
            if event.inst.subarray.tels[tel_id].optics.tel_subtype == "" and event.inst.subarray.tels[tel_id].optics.tel_type == "MST" and event.inst.subarray.tels[tel_id].optics.mirror_type == "DC":
                if tel_id not in right_tel:
                    right_tel.append(tel_id)
    return right_tel


def calibration(event, r1, dl0, dl1):
    r1.calibrate(event)
    dl0.reduce(event)
    dl1_calibrator.calibrate(event)


def set_hillas(event_info, hillas):
    arten = ["size", "cen_x", "cen_y", "length", "width", "r", "phi", "psi", "miss", "skewness", "kurtosis"]
    #              ,       m,       m,        m,       m,   m,   rad,   rad,      m
    for i in range(len(hillas)):
        art = arten[i]
        if art not in event_info:
            event_info[art] = []
        try:
            event_info[art].append(hillas[i].value)
        except:
            event_info[art].append(hillas[i])
    return event_info


def set_mc(event_info, event, tel_id):
    arten = ["mc_E", "mc_altitude", "mc_azimuth", "mc_core_x", "mc_core_y", "mc_h_first_int", "mc_azimuth_raw", "mc_altitude_raw", "mc_azimuth_cor", "mc_altitude_cor", "mc_time_slice", "mc_refstep", "tel_id", "mc_gamma_proton"]
    #           Tev,           rad,          rad,           m,           m,                m,                 ,                  ,                 ,                  ,                ,             ,         ,                  ]
    for i in arten:
        if i not in event_info:
            event_info[i] = []
    event_info["mc_E"].append(event.mc.energy.value)
    event_info["mc_altitude"].append(event.mc.alt.value)
    event_info["mc_azimuth"].append(event.mc.az.value)
    event_info["mc_core_x"].append(event.mc.core_x.value)
    event_info["mc_core_y"].append(event.mc.core_y.value)
    event_info["mc_h_first_int"].append(event.mc.h_first_int.value)
    event_info["mc_azimuth_raw"].append(event.mc.tel[tel_id].azimuth_raw)
    event_info["mc_altitude_raw"].append(event.mc.tel[tel_id].altitude_raw)
    event_info["mc_azimuth_cor"].append(event.mc.tel[tel_id].azimuth_cor)
    event_info["mc_altitude_cor"].append(event.mc.tel[tel_id].altitude_cor)
    event_info["mc_time_slice"].append(event.mc.tel[tel_id].time_slice)
    event_info["mc_refstep"].append(event.mc.tel[tel_id].meta['refstep'])
    # event_info["camera_rotation_angle"].append(event.inst.camera_rotation_angle[tel_id])
    event_info["tel_id"].append(tel_id)
    event_info["mc_gamma_proton"].append(event.mc.shower_primary_id)
    return event_info


def set_mc_header(info):
    mc_header = {}
    arten = ["mc_spectral_index", "mc_obsheight", "mc_num_showers", "mc_num_use", "mc_core_pos_mode", "mc_core_range_X", "mc_core_range_Y", "mc_alt_range_Min", "mc_alt_range_Max", "mc_az_range_Min", "mc_az_range_Max", "mc_viewcone_Min", "mc_viewcone_Max", "mc_E_range_Min", "mc_E_range_Max", "mc_diffuse", "mc_injection_height"]
    for i in arten:
        if i not in mc_header:
            mc_header[i] = info[i][0]
    return mc_header


def set_tp(save_info, cleaning_mask, mc):
    arten = ["Reinheit", "Effizienz", "Genauigkeit", "TP", "FP", "TN", "FN"]
    for i in arten:
        if i not in save_info:
            save_info[i] = []
    tp = 0
    fp = 0
    tn = 0
    fn = 0
    for i in range(len(cleaning_mask)):
        if cleaning_mask[i] == False:
            if mc[i] == 0:
                tn += 1
            else:
                fn += 1
        else:
            if mc[i] == 0:
                fp += 1
            else:
                tp += 1
    save_info["Reinheit"].append(tp / (tp + fp))
    save_info["Effizienz"].append(tp / (tp + fn))
    save_info["Genauigkeit"].append((tp + tn) / (tp + tn + fp + fn))

    save_info["TP"].append(tp)
    save_info["FP"].append(fp)
    save_info["TN"].append(tn)
    save_info["FN"].append(fn)
    return save_info


def process_event(event, dl1, infos_save):
    Keylist = ["size", "cen_x", "cen_y", "length", "width", "r", "phi", "psi", "miss", "skewness", "kurtosis", "Reinheit", "Effizienz", "Genauigkeit"]
    calibration(event, r1, dl0, dl1)

    eins_drinne = {}

    err = False
    image = event.dl1.tel[tel_id].image[0]
    geom = event.inst.subarray.tel[tel_id].camera

    for i in infos_save.keys():
        for j in infos_save[i].keys():
            if i not in eins_drinne.keys():
                eins_drinne[i] = {}
            if j not in eins_drinne[i].keys():
                eins_drinne[i][j] = False
            if i == "o":
                event_info = {}
                event_info = set_mc(event_info, event, tel_id)
                infos_save[i][j].append(event_info)
                eins_drinne[i][j] = True
                continue

            cleaning_mask = tailcuts_clean(geom, image, picture_thresh=i, boundary_thresh=j)
            if len(image[cleaning_mask]) == 0:
                continue
            clean = image.copy()
            clean[~cleaning_mask] = 0.0
            event_info = {}
            hillas = hillas_parameters(geom, image=clean)
            event_info = set_hillas(event_info, hillas)
            event_info = set_mc(event_info, event, tel_id)
            event_info = set_tp(event_info, cleaning_mask, event.mc.tel[tel_id].photo_electron_image)
            for key in Keylist:
                if key not in event_info:
                    err = True
                else:
                    drinne = False
                    wert = event_info[key][0]
                    if wert >= 0:
                        drinne = True
                    elif wert < 0:
                        drinne = True
                    if drinne is False:
                        err = True
            if err:
                continue
            infos_save[i][j].append(event_info)
            eins_drinne[i][j] = True

    return infos_save


# Set Filename and number of File
Filename = "../Master_Daten/PROD3/LaPalma/gamma/gamma_20deg_0deg_run2___cta-prod3-lapalma3-2147m-LaPalma_szip_0.hdf5"
nummer = str(0)
if len(sys.argv) == 3:
    nummer = sys.argv[1]
    Filename = sys.argv[2]

# Set Calibrator
r1 = HessioR1Calibrator(None, None)
dl0 = CameraDL0Reducer(None, None)
dl1_calibrator = CameraDL1Calibrator(
    config=None,
    tool=None,
    extractor=None,
)

infos_save = {}
infos_save["o"] = {}
infos_save["o"]["o"] = []
infos_save[7] = {}
infos_save[7][2] = []

os.system('mkdir -p preprocess_pickle')

f = h5py.File(Filename, 'r+')
f_e = f.get('events')

mc_header = {}
E_Event = []

anzahl = 0
pbar = tqdm(total=len(f_e["mc_E"]))
for i in range(len(f_e["mc_E"])):
    E_Event.append(f_e["mc_E"][i])
    if anzahl == 0:
        mc_header = set_mc_header(f.get('mc_header'))
        print(mc_header)

    #infos_save = process_event(event, dl1_calibrator, infos_save)

    # print Process time
    anzahl += 1
    if anzahl % 25 == 0:
        pbar.update(25)
pbar.update(anzahl % 25)
pbar.close()

# save Data
for i in infos_save.keys():
    for j in infos_save[i].keys():
        data_container = {"info": infos_save[i][j], "E_events": E_Event, "E_cut": E_info[i][j], "mc_header": mc_header}
        pickle.dump(data_container, open("preprocess_pickle/F" + str(nummer) + "_PT" + str(i) + "_BT" + str(j) + "_ergebnisse.pickle", "wb"))