diff --git a/data/download_data_grand.py b/data/download_data_grand.py index 3bb1e5f0..2509ef3d 100644 --- a/data/download_data_grand.py +++ b/data/download_data_grand.py @@ -18,12 +18,25 @@ #LINK_MODEL = "https://forge.in2p3.fr/attachments/download/133380/grand_model_2207.tar.gz" #FILE_MODEL = "grand_model_2207.tar.gz" -#LINK_MODEL = "https://forge.in2p3.fr/attachments/download/201909/grand_model_2306.tar.gz" -#LINK_MODEL = "https://forge.in2p3.fr/attachments/download/251637/grand_model_190224.tar.gz" -LINK_MODEL = "https://forge.in2p3.fr/attachments/download/340302/grand_model_141124.tar.gz" +LINK_MODEL = "https://forge.in2p3.fr/attachments/download/211902/grand_model_2306.tar.gz" FILE_MODEL = LINK_MODEL.split("/")[-1] +# class MyProgressBar(): +# def __init__(self): +# self.pbar = None +# +# def __call__(self, block_num, block_size, total_size): +# if not self.pbar: +# self.pbar=progressbar.ProgressBar(maxval=total_size) +# self.pbar.start() +# +# downloaded = block_num * block_size +# if downloaded < total_size: +# self.pbar.update(downloaded) +# else: +# self.pbar.finish() +# 1- test if download is necessary if os.path.exists(grand_add_path_data('detector')): print("==============================") print('Skip download data model') @@ -33,7 +46,7 @@ # 2- download print("==============================") -print("Download data model (~ 1GB) for GRAND, please wait ...") +print("Download data model (~ 452 MB) for GRAND, please wait ...") try: request.urlretrieve(LINK_MODEL, tar_file) print("Successfully downloaded") @@ -53,4 +66,4 @@ print(f"Extract failed '{tar_file}'") sys.exit(1) print("data model available in grand/data directory !") -sys.exit(0) \ No newline at end of file +sys.exit(0) \ No newline at end of file diff --git a/data/download_new_RFchain.py b/data/download_new_RFchain.py new file mode 100755 index 00000000..80099a14 --- /dev/null +++ b/data/download_new_RFchain.py @@ -0,0 +1,72 @@ +#! /usr/bin/env python3 + +''' +Created on 19 juil. 2022 + +@author: Jean-Marc Colley, CNRS/IN2P3/LPNHE + +''' +import tarfile +import os +import sys +import os.path as osp +from urllib import request +#TODO: add progressbar to grand lib +#import progressbar + +from grand import GRAND_DATA_PATH, grand_add_path_data + +#LINK_MODEL = "https://forge.in2p3.fr/attachments/download/133380/grand_model_2207.tar.gz" +#FILE_MODEL = "grand_model_2207.tar.gz" +#LINK_MODEL = "https://forge.in2p3.fr/attachments/download/201909/grand_model_2306.tar.gz" +#LINK_MODEL = "https://forge.in2p3.fr/attachments/download/251637/grand_model_190224.tar.gz" +LINK_MODEL = "https://forge.in2p3.fr/attachments/download/312643/RFchain_V3.tar.gz" +FILE_MODEL = LINK_MODEL.split("/")[-1] +# class MyProgressBar(): +# def __init__(self): +# self.pbar = None +# +# def __call__(self, block_num, block_size, total_size): +# if not self.pbar: +# self.pbar=progressbar.ProgressBar(maxval=total_size) +# self.pbar.start() +# +# downloaded = block_num * block_size +# if downloaded < total_size: +# self.pbar.update(downloaded) +# else: +# self.pbar.finish() + + +# 1- test if download is necessary +if os.path.exists(grand_add_path_data('detector/RFchain_v2')): + print("==============================") + print('Skip download data model') + sys.exit(0) + +GRAND_DATA_PATH_1 = osp.join(grand_add_path_data('detector')) +tar_file = osp.join(GRAND_DATA_PATH_1, FILE_MODEL) + +# 2- download +print("==============================") +print("Download new RFchain model (~ 450 KB) for GRAND, please wait ...") +try: + request.urlretrieve(LINK_MODEL, tar_file) + print("Successfully downloaded") +except: + print(f"download failed {LINK_MODEL}") + sys.exit(1) + +# 3- extract +print("==============================") +print('Extract tar file') +try: + my_tar = tarfile.open(tar_file) + my_tar.extractall(grand_add_path_data('detector')) + my_tar.close() + os.remove(tar_file) # delete zipped file are extraction. +except: + print(f"Extract failed '{tar_file}'") + sys.exit(1) +print("data model available in grand/data/detector directory !") +sys.exit(0) diff --git a/env/setup.sh b/env/setup.sh index 4503d654..d4ed77c7 100755 --- a/env/setup.sh +++ b/env/setup.sh @@ -10,5 +10,6 @@ cd $script_full_path # download data model for GRAND data/download_data_grand.py +data/download_new_RFchain.py cd $call_path \ No newline at end of file diff --git a/grand/sim/detector/rf_chain.py b/grand/sim/detector/rf_chain.py index 87de8558..8e98bbca 100755 --- a/grand/sim/detector/rf_chain.py +++ b/grand/sim/detector/rf_chain.py @@ -149,7 +149,7 @@ def set_out_freq_mhz(self, freqs_mhz): self.size_sig = (self.nb_freqs - 1) * 2 class MatchingNetwork(GenericProcessingDU): - + def __init__(self): """ @@ -188,7 +188,7 @@ def _set_name_data_file(self, axis): return grand_add_path_data(filename) def compute_for_freqs(self, freqs_mhz): - + logger.debug(f"{self.sparams[0].shape}") self.set_out_freq_mhz(freqs_mhz) assert self.nb_freqs > 0 @@ -263,7 +263,7 @@ def compute_for_freqs(self, freqs_mhz): class gaa_frontend0db(GenericProcessingDU): - + def __init__(self): """ @@ -302,7 +302,7 @@ def _set_name_data_file(self, axis): return grand_add_path_data(filename) def compute_for_freqs(self, freqs_mhz): - + logger.debug(f"{self.sparams[0].shape}") self.set_out_freq_mhz(freqs_mhz) assert self.nb_freqs > 0 @@ -376,7 +376,7 @@ def compute_for_freqs(self, freqs_mhz): self.ABCD_matrix[:] = ABCD_matrix # this is an A-matrix represented by [A] in the document. - + class LowNoiseAmplifier(GenericProcessingDU): """ @@ -855,7 +855,7 @@ def compute_for_freqs(self, freqs_mhz): self.ABCD_matrix[:] = s2abcd(self.s11, self.s21, self.s12, self.s22) * denorm_factor # this is an A-matrix represented by [A] in the document. ################################################################################################# - + class Rfchain_elements_db(GenericProcessingDU): def __init__(self, filename="test2.s2p"): super().__init__() @@ -892,34 +892,34 @@ def compute_for_freqs(self, freqs_mhz): dbs11 = self.sparams[:, 1] phs11 = np.deg2rad(self.sparams[:, 2]) res11, ims11 = db2reim(dbs11, phs11) - self.dbs11[:] = interpol_at_new_x(freqs_in, dbs11, self.freqs_mhz) - self.s11[:] = interpol_at_new_x(freqs_in, res11, self.freqs_mhz) - self.s11[:] += 1j * interpol_at_new_x(freqs_in, ims11, self.freqs_mhz) + self.dbs11[:] = interpol_at_new_x(freqs_in, dbs11, self.freqs_mhz) + self.s11[:] = interpol_at_new_x(freqs_in, res11, self.freqs_mhz) + self.s11[:] += 1j * interpol_at_new_x(freqs_in, ims11, self.freqs_mhz) # ----- S21 dbs21 = self.sparams[:, 3] phs21 = np.deg2rad(self.sparams[:, 4]) res21, ims21 = db2reim(dbs21, phs21) - self.dbs21[:] = interpol_at_new_x(freqs_in, dbs21, self.freqs_mhz) - self.s21[:] = interpol_at_new_x(freqs_in, res21, self.freqs_mhz) - self.s21[:] += 1j * interpol_at_new_x(freqs_in, ims21, self.freqs_mhz) + self.dbs21[:] = interpol_at_new_x(freqs_in, dbs21, self.freqs_mhz) + self.s21[:] = interpol_at_new_x(freqs_in, res21, self.freqs_mhz) + self.s21[:] += 1j * interpol_at_new_x(freqs_in, ims21, self.freqs_mhz) # ----- S12 dbs12 = self.sparams[:, 5] phs12 = np.deg2rad(self.sparams[:, 6]) res12, ims12 = db2reim(dbs12, phs12) - self.dbs12[:] = interpol_at_new_x(freqs_in, dbs12, self.freqs_mhz) - self.s12[:] = interpol_at_new_x(freqs_in, res12, self.freqs_mhz) - self.s12[:] += 1j * interpol_at_new_x(freqs_in, ims12, self.freqs_mhz) + self.dbs12[:] = interpol_at_new_x(freqs_in, dbs12, self.freqs_mhz) + self.s12[:] = interpol_at_new_x(freqs_in, res12, self.freqs_mhz) + self.s12[:] += 1j * interpol_at_new_x(freqs_in, ims12, self.freqs_mhz) # ----- S22 dbs22 = self.sparams[:][:, 7] phs22 = np.deg2rad(self.sparams[:, 8]) res22, ims22 = db2reim(dbs22, phs22) - self.dbs22[:] = interpol_at_new_x(freqs_in, dbs22, self.freqs_mhz) - self.s22[:] = interpol_at_new_x(freqs_in, res22, self.freqs_mhz) - self.s22[:] += 1j * interpol_at_new_x(freqs_in, ims22, self.freqs_mhz) + self.dbs22[:] = interpol_at_new_x(freqs_in, dbs22, self.freqs_mhz) + self.s22[:] = interpol_at_new_x(freqs_in, res22, self.freqs_mhz) + self.s22[:] += 1j * interpol_at_new_x(freqs_in, ims22, self.freqs_mhz) denorm_factor = np.array([[1, 50], [1/50., 1]]) # denormalizing factor for XYZ arms denorm_factor = denorm_factor[..., np.newaxis, np.newaxis] - # for all three ports. shape should be (2, 2, ant ports, nb_freqs) + # for all three ports. shape should be (2, 2, ant ports, nb_freqs) ABCD_matrix = s2abcd(self.s11, self.s21, self.s12, self.s22) ABCD_matrix *= denorm_factor # denormalizing factor for XYZ arms self.ABCD_matrix[:] = ABCD_matrix @@ -962,45 +962,45 @@ def compute_for_freqs(self, freqs_mhz): dbs11 = self.sparams[:, 1] phs11 = self.sparams[:, 2] res11, ims11 = db2reim(dbs11, phs11) - self.dbs11[:] = interpol_at_new_x(freqs_in, dbs11, self.freqs_mhz) - self.s11[:] = interpol_at_new_x(freqs_in, res11, self.freqs_mhz) - self.s11[:] += 1j * interpol_at_new_x(freqs_in, ims11, self.freqs_mhz) + self.dbs11[:] = interpol_at_new_x(freqs_in, dbs11, self.freqs_mhz) + self.s11[:] = interpol_at_new_x(freqs_in, res11, self.freqs_mhz) + self.s11[:] += 1j * interpol_at_new_x(freqs_in, ims11, self.freqs_mhz) # ----- S21 dbs21 = self.sparams[:, 3] phs21 = self.sparams[:, 4] res21, ims21 = db2reim(dbs21, phs21) - self.dbs21[:] = interpol_at_new_x(freqs_in, dbs21, self.freqs_mhz) - self.s21[:] = interpol_at_new_x(freqs_in, res21, self.freqs_mhz) - self.s21[:] += 1j * interpol_at_new_x(freqs_in, ims21, self.freqs_mhz) + self.dbs21[:] = interpol_at_new_x(freqs_in, dbs21, self.freqs_mhz) + self.s21[:] = interpol_at_new_x(freqs_in, res21, self.freqs_mhz) + self.s21[:] += 1j * interpol_at_new_x(freqs_in, ims21, self.freqs_mhz) # ----- S12 dbs12 = self.sparams[:, 5] phs12 = self.sparams[:, 6] res12, ims12 = db2reim(dbs12, phs12) - self.dbs12[:] = interpol_at_new_x(freqs_in, dbs12, self.freqs_mhz) - self.s12[:] = interpol_at_new_x(freqs_in, res12, self.freqs_mhz) - self.s12[:] += 1j * interpol_at_new_x(freqs_in, ims12, self.freqs_mhz) + self.dbs12[:] = interpol_at_new_x(freqs_in, dbs12, self.freqs_mhz) + self.s12[:] = interpol_at_new_x(freqs_in, res12, self.freqs_mhz) + self.s12[:] += 1j * interpol_at_new_x(freqs_in, ims12, self.freqs_mhz) # ----- S22 dbs22 = self.sparams[:][:, 7] phs22 = self.sparams[:, 8] res22, ims22 = db2reim(dbs22, phs22) - self.dbs22[:] = interpol_at_new_x(freqs_in, dbs22, self.freqs_mhz) - self.s22[:] = interpol_at_new_x(freqs_in, res22, self.freqs_mhz) - self.s22[:] += 1j * interpol_at_new_x(freqs_in, ims22, self.freqs_mhz) + self.dbs22[:] = interpol_at_new_x(freqs_in, dbs22, self.freqs_mhz) + self.s22[:] = interpol_at_new_x(freqs_in, res22, self.freqs_mhz) + self.s22[:] += 1j * interpol_at_new_x(freqs_in, ims22, self.freqs_mhz) denorm_factor = np.array([[1, 50], [1/50., 1]]) # denormalizing factor for XYZ arms denorm_factor = denorm_factor[..., np.newaxis, np.newaxis] - # for all three ports. shape should be (2, 2, ant ports, nb_freqs) + # for all three ports. shape should be (2, 2, ant ports, nb_freqs) ABCD_matrix = s2abcd(self.s11, self.s21, self.s12, self.s22) ABCD_matrix *= denorm_factor # denormalizing factor for XYZ arms self.ABCD_matrix[:] = ABCD_matrix - -###################################################################################### + +###################################################################################### class Rfchain_elements(GenericProcessingDU): def __init__(self, filename="test.s2p"): super().__init__() self.filename = filename - + self.sparams = np.loadtxt(self._set_name_data_file(), comments=['#', '!']) self.freqs_in = self.sparams[:, 0] / 1e6 # Hz to MHz # shape = (antenna_port, nb_freqs) @@ -1036,41 +1036,41 @@ def compute_for_freqs(self, freqs_mhz): angs11 = np.deg2rad(self.sparams[:, 2]) res11 = mags11 * np.cos(angs11) ims11 = mags11 * np.sin(angs11) - self.s11[:] = interpol_at_new_x(freqs_in, res11, self.freqs_mhz) - self.s11[:] += 1j * interpol_at_new_x(freqs_in, ims11, self.freqs_mhz) + self.s11[:] = interpol_at_new_x(freqs_in, res11, self.freqs_mhz) + self.s11[:] += 1j * interpol_at_new_x(freqs_in, ims11, self.freqs_mhz) # ----- S21 mags21 = self.sparams[:, 3] angs21 = np.deg2rad(self.sparams[:, 4]) res21 = mags21 * np.cos(angs21) ims21 = mags21 * np.sin(angs21) - self.s21[:] = interpol_at_new_x(freqs_in, res21, self.freqs_mhz) - self.s21[:] += 1j * interpol_at_new_x(freqs_in, ims21, self.freqs_mhz) + self.s21[:] = interpol_at_new_x(freqs_in, res21, self.freqs_mhz) + self.s21[:] += 1j * interpol_at_new_x(freqs_in, ims21, self.freqs_mhz) # ----- S12 mags12 = self.sparams[:, 5] angs12 = np.deg2rad(self.sparams[:, 6]) res12 = mags12 * np.cos(angs12) ims12 = mags12 * np.sin(angs12) - self.s12[:] = interpol_at_new_x(freqs_in, res12, self.freqs_mhz) - self.s12[:] += 1j * interpol_at_new_x(freqs_in, ims12, self.freqs_mhz) + self.s12[:] = interpol_at_new_x(freqs_in, res12, self.freqs_mhz) + self.s12[:] += 1j * interpol_at_new_x(freqs_in, ims12, self.freqs_mhz) # ----- S22 mags22 = self.sparams[:, 7] angs22 = np.deg2rad(self.sparams[:, 8]) res22 = mags22 * np.cos(angs22) ims22 = mags22 * np.sin(angs22) - self.s22[:] = interpol_at_new_x(freqs_in, res22, self.freqs_mhz) - self.s22[:] += 1j * interpol_at_new_x(freqs_in, ims22, self.freqs_mhz) + self.s22[:] = interpol_at_new_x(freqs_in, res22, self.freqs_mhz) + self.s22[:] += 1j * interpol_at_new_x(freqs_in, ims22, self.freqs_mhz) # for all three ports. shape should be (2, 2, ant ports, nb_freqs) denorm_factor = np.array([[1, 50], [1/50., 1]]) # denormalizing factor for XYZ arms denorm_factor = denorm_factor[..., np.newaxis, np.newaxis] - self.ABCD_matrix[:] = s2abcd(self.s11, self.s21, self.s12, self.s22) * denorm_factor + self.ABCD_matrix[:] = s2abcd(self.s11, self.s21, self.s12, self.s22) * denorm_factor ############################################################################################## class Rfchain_elements_rad(GenericProcessingDU): def __init__(self, filename="test.s2p"): super().__init__() self.filename = filename - + self.sparams = np.loadtxt(self._set_name_data_file(), comments=['#', '!']) self.freqs_in = self.sparams[:, 0] / 1e6 # Hz to MHz # shape = (antenna_port, nb_freqs) @@ -1106,34 +1106,34 @@ def compute_for_freqs(self, freqs_mhz): angs11 = self.sparams[:, 2] res11 = mags11 * np.cos(angs11) ims11 = mags11 * np.sin(angs11) - self.s11[:] = interpol_at_new_x(freqs_in, res11, self.freqs_mhz) - self.s11[:] += 1j * interpol_at_new_x(freqs_in, ims11, self.freqs_mhz) + self.s11[:] = interpol_at_new_x(freqs_in, res11, self.freqs_mhz) + self.s11[:] += 1j * interpol_at_new_x(freqs_in, ims11, self.freqs_mhz) # ----- S21 mags21 = self.sparams[:, 3] angs21 = self.sparams[:, 4] res21 = mags21 * np.cos(angs21) ims21 = mags21 * np.sin(angs21) - self.s21[:] = interpol_at_new_x(freqs_in, res21, self.freqs_mhz) - self.s21[:] += 1j * interpol_at_new_x(freqs_in, ims21, self.freqs_mhz) + self.s21[:] = interpol_at_new_x(freqs_in, res21, self.freqs_mhz) + self.s21[:] += 1j * interpol_at_new_x(freqs_in, ims21, self.freqs_mhz) # ----- S12 mags12 = self.sparams[:, 5] angs12 = self.sparams[:, 6] res12 = mags12 * np.cos(angs12) ims12 = mags12 * np.sin(angs12) - self.s12[:] = interpol_at_new_x(freqs_in, res12, self.freqs_mhz) - self.s12[:] += 1j * interpol_at_new_x(freqs_in, ims12, self.freqs_mhz) + self.s12[:] = interpol_at_new_x(freqs_in, res12, self.freqs_mhz) + self.s12[:] += 1j * interpol_at_new_x(freqs_in, ims12, self.freqs_mhz) # ----- S22 mags22 = self.sparams[:, 7] angs22 = self.sparams[:, 8] res22 = mags22 * np.cos(angs22) ims22 = mags22 * np.sin(angs22) - self.s22[:] = interpol_at_new_x(freqs_in, res22, self.freqs_mhz) - self.s22[:] += 1j * interpol_at_new_x(freqs_in, ims22, self.freqs_mhz) + self.s22[:] = interpol_at_new_x(freqs_in, res22, self.freqs_mhz) + self.s22[:] += 1j * interpol_at_new_x(freqs_in, ims22, self.freqs_mhz) # for all three ports. shape should be (2, 2, ant ports, nb_freqs) denorm_factor = np.array([[1, 50], [1/50., 1]]) # denormalizing factor for XYZ arms denorm_factor = denorm_factor[..., np.newaxis, np.newaxis] - self.ABCD_matrix[:] = s2abcd(self.s11, self.s21, self.s12, self.s22) * denorm_factor + self.ABCD_matrix[:] = s2abcd(self.s11, self.s21, self.s12, self.s22) * denorm_factor ########################################################################################### @@ -1163,8 +1163,8 @@ def compute_for_freqs(self, freqs_mhz): dbs = self.sparams[:, 1] phs = np.deg2rad(self.sparams[:, 2]) res, ims = db2reim(dbs, phs) - self.s[:] = interpol_at_new_x(freqs_in, res, self.freqs_mhz) - self.s[:] += 1j * interpol_at_new_x(freqs_in, ims, self.freqs_mhz) + self.s[:] = interpol_at_new_x(freqs_in, res, self.freqs_mhz) + self.s[:] += 1j * interpol_at_new_x(freqs_in, ims, self.freqs_mhz) # Calculation of Zload (Zload = balun+200ohm + ADchip) self.Z_load[:] = 50 * (1 + self.s) / (1 - self.s) @@ -1266,11 +1266,11 @@ def compute_for_freqs(self, freqs_mhz): assert self.lna.nb_freqs > 0 assert self.lna.ABCD_matrix.shape[-1] > 0 assert self.lna.nb_freqs==self.balun1.nb_freqs - + assert self.matcnet.nb_freqs > 0 assert self.matcnet.ABCD_matrix.shape[-1] > 0 assert self.matcnet.nb_freqs==self.balun1.nb_freqs - + self.Z_ant = np.zeros((3, self.nb_freqs), dtype=np.complex64) self.Z_in = np.zeros((3, self.nb_freqs), dtype=np.complex64) self.V_out_RFchain = np.zeros((3, self.nb_freqs), dtype=np.complex64) @@ -1281,20 +1281,20 @@ def compute_for_freqs(self, freqs_mhz): # Note that this is a matrix multiplication # Associative property of matrix multiplication is used, ie. (AB)C = A(BC) # Make sure to multiply in this order: balun1 * matching_network * lna.ABCD_matrix * cable.ABCD_matrix * vgaf.ABCD_matrix - + MM1 = matmul(self.balun1.ABCD_matrix, self.matcnet.ABCD_matrix) MM2 = matmul(MM1, self.lna.ABCD_matrix) MM3 = matmul(self.cable.ABCD_matrix, self.vgaf.ABCD_matrix) self.total_ABCD_matrix[:] = matmul(MM2, MM3) - + # Calculation of Z_in (this is the total impedence of the RF chain excluding antenna arm. see page 50 of the document.) self.Z_load = self.zload.Z_load[np.newaxis, :] # shape (nfreq) --> (1,nfreq) to broadcast with components of ABCD_matrix with shape (2,2,ports,nfreq). self.Z_in[:] = (self.total_ABCD_matrix[0,0] * self.Z_load + self.total_ABCD_matrix[0,1])/(self.total_ABCD_matrix[1,0] * self.Z_load + self.total_ABCD_matrix[1,1]) #self.Z_in[:] = (self.total_ABCD_matrix[0,0] * self.Z_load + self.total_ABCD_matrix[0,1])/(self.total_ABCD_matrix[1,0] * self.Z_load + self.total_ABCD_matrix[1,1]) - + # Once Z_in is calculated, calculate the final total_ABCD_matrix including Balun2. self.total_ABCD_matrix[:] = matmul(self.total_ABCD_matrix, self.balun2.ABCD_matrix) - #self.total_ABCD_matrix[:] = matmul(self.total_ABCD_matrix, self.balun2.ABCD_matrix) + #self.total_ABCD_matrix[:] = matmul(self.total_ABCD_matrix, self.balun2.ABCD_matrix) # Antenna Impedance. filename = os.path.join("detector", "RFchain_v2", "Z_ant_3.2m.csv") @@ -1313,7 +1313,7 @@ def vout_f(self, voc_f): Input: Voc_f (in frequency domain) Output: Voltage after RF chain in frequency domain. Make sure to run self.compute_for_freqs() before calling this method. - RK Note: name 'vout_f' is a placeholder. Change it with something better. + RK Note: name 'vout_f' is a placeholder. Change it with something better. """ assert voc_f.shape==self.Z_in.shape # shape = (nports, nfreqs) @@ -1393,21 +1393,21 @@ def compute_for_freqs(self, freqs_mhz): # Note that this is a matrix multiplication # Associative property of matrix multiplication is used, ie. (AB)C = A(BC) # Make sure to multiply in this order: lna.ABCD_matrix * balun1.ABCD_matrix * cable.ABCD_matrix * vgaf.ABCD_matrix - + MM1 = matmul(self.balun1.ABCD_matrix, self.matcnet.ABCD_matrix) MM2 = matmul(MM1, self.lna.ABCD_matrix) MM3 = matmul(self.cable.ABCD_matrix, self.vgaf.ABCD_matrix) self.total_ABCD_matrix[:] = matmul(MM2, MM3) - + MMM1 = matmul(self.balun1.ABCD_matrix, self.matcnet.ABCD_matrix) self.total_ABCD_matrix_nut[:] = matmul(MMM1, self.lna.ABCD_matrix) - + # Calculation of Z_in (this is the total impedence of the RF chain excluding antenna arm. see page 50 of the document.) self.Z_load = self.zload.Z_load[np.newaxis, :] # shape (nfreq) --> (1,nfreq) to broadcast with components of ABCD_matrix with shape (2,2,ports,nfreq). self.Z_in[:] = (self.total_ABCD_matrix[0,0] * self.Z_load + self.total_ABCD_matrix[0,1])/(self.total_ABCD_matrix[1,0] * self.Z_load + self.total_ABCD_matrix[1,1]) # Once Z_in is calculated, calculate the final total_ABCD_matrix including Balun2. - #self.total_ABCD_matrix[:] = matmul(self.total_ABCD_matrix, self.balun2.ABCD_matrix) + #self.total_ABCD_matrix[:] = matmul(self.total_ABCD_matrix, self.balun2.ABCD_matrix) # Antenna Impedance. filename = os.path.join("detector", "RFchain_v2", "Z_ant_3.2m.csv") @@ -1426,7 +1426,7 @@ def vout_f(self, voc_f): Input: Voc_f (in frequency domain) Output: Voltage after RF chain in frequency domain. Make sure to run self.compute_for_freqs() before calling this method. - RK Note: name 'vout_f' is a placeholder. Change it with something better. + RK Note: name 'vout_f' is a placeholder. Change it with something better. """ assert voc_f.shape==self.Z_in.shape # shape = (nports, nfreqs) @@ -1455,7 +1455,7 @@ def get_tf(self): self._total_tf = self.vout_f(np.ones((3, self.nb_freqs))) return self._total_tf -################################################################################## +################################################################################## class RFChain_gaa(GenericProcessingDU): """ @@ -1492,15 +1492,15 @@ def compute_for_freqs(self, freqs_mhz): self.V_out_RFchain = np.zeros((3, self.nb_freqs), dtype=np.complex64) self.I_out_RFchain = np.zeros((3, self.nb_freqs), dtype=np.complex64) self.total_ABCD_matrix = np.zeros(self.gaa.ABCD_matrix.shape, dtype=np.complex64) - + self.total_ABCD_matrix[:] = self.gaa.ABCD_matrix - + # Calculation of Z_in (this is the total impedence of the RF chain excluding antenna arm. see page 50 of the document.) self.Z_load = self.zload.Z_load[np.newaxis, :] # shape (nfreq) --> (1,nfreq) to broadcast with components of ABCD_matrix with shape (2,2,ports,nfreq). self.Z_in[:] = (self.total_ABCD_matrix[0,0] * self.Z_load + self.total_ABCD_matrix[0,1])/(self.total_ABCD_matrix[1,0] * self.Z_load + self.total_ABCD_matrix[1,1]) # Once Z_in is calculated, calculate the final total_ABCD_matrix including Balun2. - #self.total_ABCD_matrix[:] = matmul(self.total_ABCD_matrix, self.balun2.ABCD_matrix) + #self.total_ABCD_matrix[:] = matmul(self.total_ABCD_matrix, self.balun2.ABCD_matrix) # Antenna Impedance. filename = os.path.join("detector", "RFchain_v2", "Z_ant_3.2m.csv") @@ -1519,7 +1519,7 @@ def vout_f(self, voc_f): Input: Voc_f (in frequency domain) Output: Voltage after RF chain in frequency domain. Make sure to run self.compute_for_freqs() before calling this method. - RK Note: name 'vout_f' is a placeholder. Change it with something better. + RK Note: name 'vout_f' is a placeholder. Change it with something better. """ assert voc_f.shape==self.Z_in.shape # shape = (nports, nfreqs) @@ -1532,7 +1532,7 @@ def vout_f(self, voc_f): ABCD_matrix_2port = self.total_ABCD_matrix[:,:,i,:] ABCD_matrix_1port = np.moveaxis(ABCD_matrix_1port, -1, 0) # (2,2,nfreqs) --> (nfreqs,2,2), to compute inverse of ABCD_matrix using np.linalg.inv. ABCD_matrix_1port_inv = np.linalg.inv(ABCD_matrix_1port) - + self.V_out_RFchain[i] = 2/(self.total_ABCD_matrix[0,0,i,:] + self.total_ABCD_matrix[0,1,i,:]/50 + self.total_ABCD_matrix[1,0,i,:]*50 + self.total_ABCD_matrix[1,1,i,:]) #V_out_RFchain = ABCD_matrix_1port_inv[:,0,0]*self.V_in_balunA[i] + ABCD_matrix_1port_inv[:,0,1]*self.I_in_balunA[i] I_out_RFchain = ABCD_matrix_1port_inv[:,1,0]*self.V_in_balunA[i] + ABCD_matrix_1port_inv[:,1,1]*self.I_in_balunA[i] @@ -1554,7 +1554,7 @@ def get_tf(self): return self._total_tf ########################################################################################## ########################################################################################### - + class RFChain_Balun1(GenericProcessingDU): """ Facade for all elements in RF chain @@ -1606,20 +1606,20 @@ def compute_for_freqs(self, freqs_mhz): # Note that this is a matrix multiplication # Associative property of matrix multiplication is used, ie. (AB)C = A(BC) # Make sure to multiply in this order: lna.ABCD_matrix * balun1.ABCD_matrix * cable.ABCD_matrix * vgaf.ABCD_matrix - + MM1 = matmul(self.balun1.ABCD_matrix, self.matcnet.ABCD_matrix) MM2 = matmul(MM1, self.lna.ABCD_matrix) MM3 = matmul(self.cable.ABCD_matrix, self.vgaf.ABCD_matrix) self.total_ABCD_matrix[:] = matmul(MM2, MM3) - + self.total_ABCD_matrix_nut[:] = self.balun1.ABCD_matrix - + # Calculation of Z_in (this is the total impedence of the RF chain excluding antenna arm. see page 50 of the document.) self.Z_load = self.zload.Z_load[np.newaxis, :] # shape (nfreq) --> (1,nfreq) to broadcast with components of ABCD_matrix with shape (2,2,ports,nfreq). self.Z_in[:] = (self.total_ABCD_matrix[0,0] * self.Z_load + self.total_ABCD_matrix[0,1])/(self.total_ABCD_matrix[1,0] * self.Z_load + self.total_ABCD_matrix[1,1]) # Once Z_in is calculated, calculate the final total_ABCD_matrix including Balun2. - #self.total_ABCD_matrix[:] = matmul(self.total_ABCD_matrix, self.balun2.ABCD_matrix) + #self.total_ABCD_matrix[:] = matmul(self.total_ABCD_matrix, self.balun2.ABCD_matrix) # Antenna Impedance. filename = os.path.join("detector", "RFchain_v2", "Z_ant_3.2m.csv") @@ -1638,7 +1638,7 @@ def vout_f(self, voc_f): Input: Voc_f (in frequency domain) Output: Voltage after RF chain in frequency domain. Make sure to run self.compute_for_freqs() before calling this method. - RK Note: name 'vout_f' is a placeholder. Change it with something better. + RK Note: name 'vout_f' is a placeholder. Change it with something better. """ assert voc_f.shape==self.Z_in.shape # shape = (nports, nfreqs) @@ -1667,9 +1667,9 @@ def get_tf(self): self._total_tf = self.vout_f(np.ones((3, self.nb_freqs))) return self._total_tf -################################################################################## +################################################################################## ########################################################################################### - + class RFChain_Match_net(GenericProcessingDU): """ Facade for all elements in RF chain @@ -1721,20 +1721,20 @@ def compute_for_freqs(self, freqs_mhz): # Note that this is a matrix multiplication # Associative property of matrix multiplication is used, ie. (AB)C = A(BC) # Make sure to multiply in this order: lna.ABCD_matrix * balun1.ABCD_matrix * cable.ABCD_matrix * vgaf.ABCD_matrix - + MM1 = matmul(self.balun1.ABCD_matrix, self.matcnet.ABCD_matrix) MM2 = matmul(MM1, self.lna.ABCD_matrix) MM3 = matmul(self.cable.ABCD_matrix, self.vgaf.ABCD_matrix) self.total_ABCD_matrix[:] = matmul(MM2, MM3) - + self.total_ABCD_matrix_nut[:] = matmul(self.balun1.ABCD_matrix, self.matcnet.ABCD_matrix) - + # Calculation of Z_in (this is the total impedence of the RF chain excluding antenna arm. see page 50 of the document.) self.Z_load = self.zload.Z_load[np.newaxis, :] # shape (nfreq) --> (1,nfreq) to broadcast with components of ABCD_matrix with shape (2,2,ports,nfreq). self.Z_in[:] = (self.total_ABCD_matrix[0,0] * self.Z_load + self.total_ABCD_matrix[0,1])/(self.total_ABCD_matrix[1,0] * self.Z_load + self.total_ABCD_matrix[1,1]) # Once Z_in is calculated, calculate the final total_ABCD_matrix including Balun2. - #self.total_ABCD_matrix[:] = matmul(self.total_ABCD_matrix, self.balun2.ABCD_matrix) + #self.total_ABCD_matrix[:] = matmul(self.total_ABCD_matrix, self.balun2.ABCD_matrix) # Antenna Impedance. filename = os.path.join("detector", "RFchain_v2", "Z_ant_3.2m.csv") @@ -1753,7 +1753,7 @@ def vout_f(self, voc_f): Input: Voc_f (in frequency domain) Output: Voltage after RF chain in frequency domain. Make sure to run self.compute_for_freqs() before calling this method. - RK Note: name 'vout_f' is a placeholder. Change it with something better. + RK Note: name 'vout_f' is a placeholder. Change it with something better. """ assert voc_f.shape==self.Z_in.shape # shape = (nports, nfreqs) @@ -1782,7 +1782,7 @@ def get_tf(self): self._total_tf = self.vout_f(np.ones((3, self.nb_freqs))) return self._total_tf -################################################################################## +################################################################################## class RFChain_Cable_Connectors(GenericProcessingDU): """ @@ -1835,22 +1835,22 @@ def compute_for_freqs(self, freqs_mhz): # Note that this is a matrix multiplication # Associative property of matrix multiplication is used, ie. (AB)C = A(BC) # Make sure to multiply in this order: lna.ABCD_matrix * balun1.ABCD_matrix * cable.ABCD_matrix * vgaf.ABCD_matrix - + MM1 = matmul(self.balun1.ABCD_matrix, self.matcnet.ABCD_matrix) MM2 = matmul(MM1, self.lna.ABCD_matrix) MM3 = matmul(self.cable.ABCD_matrix, self.vgaf.ABCD_matrix) self.total_ABCD_matrix[:] = matmul(MM2, MM3) - + MMM1 = matmul(self.balun1.ABCD_matrix, self.matcnet.ABCD_matrix) MMM2 = matmul(MMM1, self.lna.ABCD_matrix) self.total_ABCD_matrix_nut[:] = matmul(MMM2, self.cable.ABCD_matrix) - + # Calculation of Z_in (this is the total impedence of the RF chain excluding antenna arm. see page 50 of the document.) self.Z_load = self.zload.Z_load[np.newaxis, :] # shape (nfreq) --> (1,nfreq) to broadcast with components of ABCD_matrix with shape (2,2,ports,nfreq). self.Z_in[:] = (self.total_ABCD_matrix[0,0] * self.Z_load + self.total_ABCD_matrix[0,1])/(self.total_ABCD_matrix[1,0] * self.Z_load + self.total_ABCD_matrix[1,1]) # Once Z_in is calculated, calculate the final total_ABCD_matrix including Balun2. - #self.total_ABCD_matrix[:] = matmul(self.total_ABCD_matrix, self.balun2.ABCD_matrix) + #self.total_ABCD_matrix[:] = matmul(self.total_ABCD_matrix, self.balun2.ABCD_matrix) # Antenna Impedance. filename = os.path.join("detector", "RFchain_v2", "Z_ant_3.2m.csv") @@ -1869,7 +1869,7 @@ def vout_f(self, voc_f): Input: Voc_f (in frequency domain) Output: Voltage after RF chain in frequency domain. Make sure to run self.compute_for_freqs() before calling this method. - RK Note: name 'vout_f' is a placeholder. Change it with something better. + RK Note: name 'vout_f' is a placeholder. Change it with something better. """ assert voc_f.shape==self.Z_in.shape # shape = (nports, nfreqs) @@ -1898,7 +1898,7 @@ def get_tf(self): self._total_tf = self.vout_f(np.ones((3, self.nb_freqs))) return self._total_tf -################################################################################## +################################################################################## class RFChain_VGA(GenericProcessingDU): """ @@ -1951,21 +1951,21 @@ def compute_for_freqs(self, freqs_mhz): # Note that this is a matrix multiplication # Associative property of matrix multiplication is used, ie. (AB)C = A(BC) # Make sure to multiply in this order: lna.ABCD_matrix * balun1.ABCD_matrix * cable.ABCD_matrix * vgaf.ABCD_matrix - + MM1 = matmul(self.balun1.ABCD_matrix, self.matcnet.ABCD_matrix) MM2 = matmul(MM1, self.lna.ABCD_matrix) MM3 = matmul(self.cable.ABCD_matrix, self.vgaf.ABCD_matrix) self.total_ABCD_matrix[:] = matmul(MM2, MM3) - + # Calculation of Z_in (this is the total impedence of the RF chain excluding antenna arm. see page 50 of the document.) self.Z_load = self.zload.Z_load[np.newaxis, :] # shape (nfreq) --> (1,nfreq) to broadcast with components of ABCD_matrix with shape (2,2,ports,nfreq). self.Z_in[:] = (self.total_ABCD_matrix[0,0] * self.Z_load + self.total_ABCD_matrix[0,1])/(self.total_ABCD_matrix[1,0] * self.Z_load + self.total_ABCD_matrix[1,1]) # Once Z_in is calculated, calculate the final total_ABCD_matrix including Balun2. - #self.total_ABCD_matrix[:] = matmul(self.total_ABCD_matrix, self.balun2.ABCD_matrix) + #self.total_ABCD_matrix[:] = matmul(self.total_ABCD_matrix, self.balun2.ABCD_matrix) # Antenna Impedance. - filename = os.path.join("detector", "RFchain_v2", "Z_ant_3.2m.csv") + filename = os.path.join("detector", "RFchain_v1", "Z_ant_3.2m.csv") filename = grand_add_path_data(filename) Zant_dat = np.loadtxt(filename, delimiter=",", comments=['#', '!'], skiprows=1) freqs_in = Zant_dat[:,0] # MHz @@ -2010,7 +2010,7 @@ def get_tf(self): self._total_tf = self.vout_f(np.ones((3, self.nb_freqs))) return self._total_tf -################################################################################## +################################################################################## class RFChain_in_Balun1(GenericProcessingDU): """ Facade for all elements in RF chain @@ -2062,20 +2062,20 @@ def compute_for_freqs(self, freqs_mhz): # Note that this is a matrix multiplication # Associative property of matrix multiplication is used, ie. (AB)C = A(BC) # Make sure to multiply in this order: lna.ABCD_matrix * balun1.ABCD_matrix * cable.ABCD_matrix * vgaf.ABCD_matrix - + MM1 = matmul(self.balun1.ABCD_matrix, self.matcnet.ABCD_matrix) MM2 = matmul(MM1, self.lna.ABCD_matrix) MM3 = matmul(self.cable.ABCD_matrix, self.vgaf.ABCD_matrix) self.total_ABCD_matrix[:] = matmul(MM2, MM3) - + self.total_ABCD_matrix_nut[:] = np.ones(self.lna.ABCD_matrix.shape, dtype=np.complex64) - + # Calculation of Z_in (this is the total impedence of the RF chain excluding antenna arm. see page 50 of the document.) self.Z_load = self.zload.Z_load[np.newaxis, :] # shape (nfreq) --> (1,nfreq) to broadcast with components of ABCD_matrix with shape (2,2,ports,nfreq). self.Z_in[:] = (self.total_ABCD_matrix[0,0] * self.Z_load + self.total_ABCD_matrix[0,1])/(self.total_ABCD_matrix[1,0] * self.Z_load + self.total_ABCD_matrix[1,1]) # Once Z_in is calculated, calculate the final total_ABCD_matrix including Balun2. - #self.total_ABCD_matrix[:] = matmul(self.total_ABCD_matrix, self.balun2.ABCD_matrix) + #self.total_ABCD_matrix[:] = matmul(self.total_ABCD_matrix, self.balun2.ABCD_matrix) # Antenna Impedance. filename = os.path.join("detector", "RFchain_v2", "Z_ant_3.2m.csv") @@ -2094,7 +2094,7 @@ def vout_f(self, voc_f): Input: Voc_f (in frequency domain) Output: Voltage after RF chain in frequency domain. Make sure to run self.compute_for_freqs() before calling this method. - RK Note: name 'vout_f' is a placeholder. Change it with something better. + RK Note: name 'vout_f' is a placeholder. Change it with something better. """ assert voc_f.shape==self.Z_in.shape # shape = (nports, nfreqs) @@ -2125,4 +2125,4 @@ def get_tf(self): self._total_tf = self.vout_f(np.ones((3, self.nb_freqs))) return self._total_tf -################################################################################## \ No newline at end of file +################################################################################## \ No newline at end of file diff --git a/sim2root/CoREASRawRoot/CoreasToRawROOT.py b/sim2root/CoREASRawRoot/CoreasToRawROOT.py index 28188f02..5cf69231 100755 --- a/sim2root/CoREASRawRoot/CoreasToRawROOT.py +++ b/sim2root/CoREASRawRoot/CoreasToRawROOT.py @@ -3,12 +3,13 @@ ## by Jelena Köhler, @jelenakhlr import sys +import numpy as np import os import glob import time #to get the unix timestamp from CorsikaInfoFuncs import * # this is in the same dir as this file sys.path.append("../Common") -import raw_root_trees as RawTrees # this is in Common. since we're in CoREASRawRoot, this is in ../Common +import sim2root.Common.raw_root_trees as RawTrees # this is in Common. since we're in CoREASRawRoot, this is in ../Common from optparse import OptionParser """ @@ -237,6 +238,15 @@ def CoreasToRawRoot(file, simID=None): ed_sum = energy_dep[:,9] + Egamma = np.sum(ed_gamma) + Eem_ion = np.sum(ed_em_ioniz) + Eem_cut = np.sum(ed_em_cut) + + # calculate electromagnetic shower energy and convert to eV + Eem = (Egamma + Eem_ion + Eem_cut) * 1e9 + print("Electromagnetic shower energy:", Eem) + + ############################################## EnergyInNeutrinos = 1. # placeholder @@ -261,6 +271,9 @@ def CoreasToRawRoot(file, simID=None): site = read_site(inp_input) latitude, longitude, altitude = read_lat_long_alt(site) + # set altitude to simulations obslevel + altitude = CorePosition[2] + ############################################################################################################################ # Part B.I.ii: Create and fill the RAW Shower Tree ############################################################################################################################ @@ -282,6 +295,7 @@ def CoreasToRawRoot(file, simID=None): RawShower.rnd_seed = RandomSeed RawShower.energy_in_neutrinos = EnergyInNeutrinos + RawShower.energy_em = [Eem] RawShower.energy_primary = [Energy] RawShower.azimuth = azimuth RawShower.zenith = zenith @@ -410,9 +424,12 @@ def CoreasToRawRoot(file, simID=None): #****** load positions ****** # the list file contains all antenna positions for each antenna ID pathAntennaList = f"{path}/SIM{simID}.list" + # store all antenna IDs in ant_IDs + # stores the actual antenna names from the coreas -list file antenna_names = antenna_positions_dict(pathAntennaList)["name"] - antenna_IDs = antenna_positions_dict(pathAntennaList)["ID"] + # if antenna name is in the form of "ant???" or "gp_???", stores the digits after "ant" or "gp_", otherwise generic counter from 1 + antenna_IDs = antenna_positions_dict(pathAntennaList)["ID"] ############################################################################################################################ # Part B.II.ii: Create and fill the RawEfield Tree @@ -443,7 +460,7 @@ def CoreasToRawRoot(file, simID=None): # the files are setup like [timestamp, x polarization, y polarization, z polarization] efield = np.loadtxt(tracefile) - timestamp = efield[:,0] * 10**9 #convert to ns + timestamp = efield[:,0] * 10**9 # convert to ns # coreas uses cgs units, so voltage is in statvolt # efield in statvolt / cm # 1 statV * 299.792458 V/statV = 1 V diff --git a/sim2root/CoREASRawRoot/Coreas_000004.root b/sim2root/CoREASRawRoot/Coreas_000004.root new file mode 100644 index 00000000..6b057cc0 Binary files /dev/null and b/sim2root/CoREASRawRoot/Coreas_000004.root differ diff --git a/sim2root/CoREASRawRoot/CorsikaInfoFuncs.py b/sim2root/CoREASRawRoot/CorsikaInfoFuncs.py index b7c8e180..293d9a40 100755 --- a/sim2root/CoREASRawRoot/CorsikaInfoFuncs.py +++ b/sim2root/CoREASRawRoot/CorsikaInfoFuncs.py @@ -186,11 +186,15 @@ def antenna_positions_dict(pathAntennaList): generic_id_counter = 1 for name in file[:, 5]: - GP_match = search(r'gp_(\d+)', name, flags=re.IGNORECASE) # Extract digits after last underscore + GP_match = search(r'gp_(\d+)', name, flags=re.IGNORECASE) # Extract digits after "ant" + ant_match = search(r'ant(\d+)', name, flags=re.IGNORECASE) # Extract digits after "ant" if GP_match: # match GP13 antennaInfo["ID"].append(int(GP_match.group(1))) # Group 1 contains the ID + elif ant_match: # match GP300 antenna names of "ant???" + antennaInfo["ID"].append(int(ant_match.group(1))) # Group 1 contains the ID + else: # Give generic IDs to antennas with other names antennaInfo["ID"].append(int(generic_id_counter)) generic_id_counter += 1 @@ -212,8 +216,9 @@ def get_antenna_position(pathAntennaList, antenna): or None if the antenna is not found. """ file = np.genfromtxt(pathAntennaList, dtype="str") - gp300_data = np.genfromtxt("GP300.list", dtype="str") + gp300_data = np.genfromtxt(pathAntennaList, dtype="str") + # currently checks the same file, so always passes # Filter data based on antenna name filtered_data = file[file[:, 5] == antenna] # Select rows where name matches gp300_filtered = gp300_data[gp300_data[:, 5] == antenna] @@ -222,7 +227,7 @@ def get_antenna_position(pathAntennaList, antenna): # Antenna not found in either file return None - # Replace x, y, z positions with the original positions + # Replace x, y, z positions with the original positions and convert to m x = gp300_filtered[0][2].astype(float) * 10**-2 y = gp300_filtered[0][3].astype(float) * 10**-2 z = gp300_filtered[0][4].astype(float) * 10**-2 @@ -245,13 +250,17 @@ def calculate_array_shift(pathAntennaList): If the arrays have different shapes, returns None. """ sim_data = np.genfromtxt(pathAntennaList, dtype="str") - gp300_data = np.genfromtxt("GP300.list", dtype="str") + gp300_data = np.genfromtxt(pathAntennaList, dtype="str") # Check if data shapes are compatible + # currently checks the same file so always passes if sim_data.shape != gp300_data.shape: print("Sim data shape is not compatible with GP300 layout!") return None + # currently checks the same file so shift is always 0 + # current version uses variable cores so no shift necessary + # utilise this process corsika sims that have core [0,0,0] and shift the antennas instead # Calculate average shift for x and y coordinates shift_x = np.mean(gp300_data[:, 2].astype(float) * 10**-2 - sim_data[:, 2].astype(float) * 10**-2) shift_y = np.mean(gp300_data[:, 3].astype(float) * 10**-2 - sim_data[:, 3].astype(float) * 10**-2) diff --git a/sim2root/Common/IllustrateSimPipe.py b/sim2root/Common/IllustrateSimPipe.py index 435516b6..ee9c0163 100644 --- a/sim2root/Common/IllustrateSimPipe.py +++ b/sim2root/Common/IllustrateSimPipe.py @@ -180,8 +180,7 @@ def plot_traces_all_levels(directory, t_0_shift=False): trace_voltage = np.asarray(tvoltage_l0.trace, dtype=np.float32) trace_ADC_L1= np.asarray(tadc_l1.trace_ch, dtype=np.float32) trace_efield_L1= np.asarray(tefield_l1.trace, dtype=np.float32) - du_id = np.asarray(tefield_l0.du_id) # MT: used for printing info and saving in voltage tree. - # JK: du_id is currently unused - TODO + du_id = np.asarray(tefield_l0.du_id) # t0 calculations @@ -216,8 +215,8 @@ def plot_traces_all_levels(directory, t_0_shift=False): du_xyzs= np.asarray(trun_l0.du_xyz)[event_dus_indices] - # loop over all stations. - for du_idx in range(nb_du): + # loop over all stations in the event. + for du_idx in range(len(du_id)): logger.debug(f"Running DU number {du_idx}") # efield trace L0 @@ -257,11 +256,11 @@ def plot_traces_all_levels(directory, t_0_shift=False): trace_ADC_L1_time += t0_adc_L1[du_idx] trace_efield_L1_time += t0_efield_L1[du_idx] - plt.suptitle(f"event {event_number}, run {run_number}, antenna {du_idx} - Shower Time") + plt.suptitle(f"event {event_number}, run {run_number}, antenna {du_id[du_idx]} - Shower Time") savelabel = "with_t0_shift" else: logger.debug(f"NOT shifting by t0 - peaks should be at 0ns") - plt.suptitle(f"event {event_number}, run {run_number}, antenna {du_idx} - Trigger Time") + plt.suptitle(f"event {event_number}, run {run_number}, antenna {du_id[du_idx]} - Trigger Time \n position {du_xyzs[du_idx]}") savelabel = "no_t0s_shift" @@ -303,15 +302,15 @@ def plot_traces_all_levels(directory, t_0_shift=False): if t_0_shift == True: - ax1.axvline(t0_voltage_L0[du_idx], label="t0 Trigger") - ax2.axvline(t0_adc_L1[du_idx], label="t0 Trigger") - ax3.axvline(t0_efield_L0[du_idx], label="t0 Trigger") - ax4.axvline(t0_efield_L1[du_idx], label="t0 Trigger") + ax1.axvline(t0_voltage_L0[du_idx], label="t0 Trigger", alpha=0.2) + ax2.axvline(t0_adc_L1[du_idx], label="t0 Trigger", alpha=0.2) + ax3.axvline(t0_efield_L0[du_idx], label="t0 Trigger", alpha=0.2) + ax4.axvline(t0_efield_L1[du_idx], label="t0 Trigger", alpha=0.2) else: - ax1.axvline(0, label="t0 Trigger") - ax2.axvline(0, label="t0 Trigger") - ax3.axvline(0, label="t0 Trigger") - ax4.axvline(0, label="t0 Trigger") + ax1.axvline(0, label="t0 Trigger", alpha=0.2) + ax2.axvline(0, label="t0 Trigger", alpha=0.2) + ax3.axvline(0, label="t0 Trigger", alpha=0.2) + ax4.axvline(0, label="t0 Trigger", alpha=0.2) # Add common vertical line (assuming same time axis) for ax in [ax1, ax2, ax3, ax4]: @@ -319,7 +318,7 @@ def plot_traces_all_levels(directory, t_0_shift=False): plt.tight_layout() if(args.savefig): - plt.savefig(f"{plot_dir}/plots/IllustrateSimPipe_{run_number}_{event_number}_{du_idx}_{savelabel}.png") + plt.savefig(f"{plot_dir}/plots/IllustrateSimPipe_{run_number}_{event_number}_{du_id[du_idx]}_{savelabel}.png") plt.close(fig) else: plt.show() diff --git a/sim2root/Common/raw_root_trees.py b/sim2root/Common/raw_root_trees.py index 6ed83ced..0d51c01e 100644 --- a/sim2root/Common/raw_root_trees.py +++ b/sim2root/Common/raw_root_trees.py @@ -67,6 +67,9 @@ class RawShowerTree(MotherEventTree): ###X Primary energy (GeV) energy_primary: StdVectorListDesc = field(default=StdVectorListDesc("float")) + + ###X Electromagnetic energy (GeV) + energy_em: StdVectorListDesc = field(default=StdVectorListDesc("float")) ### Shower azimuth (deg, CR convention) _azimuth: np.ndarray = field(default_factory=lambda: np.zeros(1, np.float32)) diff --git a/sim2root/Common/sim2root.py b/sim2root/Common/sim2root.py index 3673ae3c..3708d6f2 100755 --- a/sim2root/Common/sim2root.py +++ b/sim2root/Common/sim2root.py @@ -14,8 +14,7 @@ import raw_root_trees as RawTrees # this is here in Common import grand.manage_log as mlg import matplotlib.pyplot as plt -# from scipy.ndimage.interpolation import shift #to shift the time trance for the trigger simulation -# from scipy.ndimage import shift #to shift the time trance for the trigger simulation +from scipy.ndimage import shift #to shift the time trance for the trigger simulation # specific logger definition for script because __mane__ is "__main__" ! logger = mlg.get_logger_for_script(__file__) @@ -147,23 +146,21 @@ def adjust_trigger(trace, CurrentT0s, TPre, TimeBinSize): DeltaT = TPre - trigger_time ShiftBins = (DeltaT / TimeBinSize).astype(int, copy=False) - # this is to assure that, if the maximum is found too late in the trace, we dont move outside of the original time window (normally, peaks are late in the time window, if you set the time window correctly). - mask = ShiftBins < -TPre / TimeBinSize + #this is to assure that, if the maximum is found too late in the trace, we dont move outside of the original time window (normally, peaks are late in the time window, if you set the time window correctly). + mask=ShiftBins < -TPre/TimeBinSize if mask.any(): logger.error("some elements needed to be shifted only up to the limt, tpre was too small") - ShiftBins[mask] = int(-TPre / TimeBinSize) + ShiftBins[mask]= int(-TPre/TimeBinSize) - # we cannot use np.roll, but roll makes re-appear the end of the trace at the begining if we roll to much - # we cannot use scipy shift, that lets you state what value to put for the places you roll, on a 3D array + #we cannot use use np.roll, but roll makes re-appear the end of the trace at the begining if we roll to much + #we cannot use scipy shift, that lets you state what value to put for the places you roll, on a 3D array - # TODO: There must be a better way to do this without the for loop, but i lost a morning to it and i dont have the time to develop it now. Search for strided_indexing_roll on the web for inspiration. - # for du_idx in range(trace.shape[0]): - # trace[du_idx] = shift(trace[du_idx], (0, ShiftBins[du_idx]), cval=0) + #TODO: There must be a better way to do this without the for loop, but i lost a morning to it and i dont have the time to develop it now. Search for strided_indexing_roll on the web for inspiration. for du_idx in range(trace.shape[0]): - trace_shift(trace[du_idx], ShiftBins[du_idx]) + trace[du_idx]=shift(trace[du_idx],(0,ShiftBins[du_idx]),cval=0) - # we get the correct t0 - T0s = CurrentT0s - ShiftBins * TimeBinSize + #we get the correct t0 + T0s=CurrentT0s-ShiftBins*TimeBinSize return T0s, trace @@ -622,6 +619,8 @@ def rawshower2grandroot(trawshower, gt): # ToDo: it should be a scalar on sim side gt.tshower.energy_primary = trawshower.energy_primary[0] + gt.tshower.energy_em = trawshower.energy_em[0] + ### Shower azimuth (deg, CR convention) gt.tshower.azimuth = trawshower.azimuth diff --git a/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/adc_5388-23832_L1_0000.root b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/adc_5388-23832_L1_0000.root new file mode 100644 index 00000000..ab04a8f6 Binary files /dev/null and b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/adc_5388-23832_L1_0000.root differ diff --git a/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/efield_5388-23832_L0_0000.root b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/efield_5388-23832_L0_0000.root new file mode 100644 index 00000000..81d0ffa2 Binary files /dev/null and b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/efield_5388-23832_L0_0000.root differ diff --git a/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/efield_5388-23832_L1_0000.root b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/efield_5388-23832_L1_0000.root new file mode 100644 index 00000000..c00023c5 Binary files /dev/null and b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/efield_5388-23832_L1_0000.root differ diff --git a/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/run_0_L0_0000.root b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/run_0_L0_0000.root new file mode 100644 index 00000000..7a56e9bb Binary files /dev/null and b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/run_0_L0_0000.root differ diff --git a/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/run_0_L1_0000.root b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/run_0_L1_0000.root new file mode 100644 index 00000000..13012f2e Binary files /dev/null and b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/run_0_L1_0000.root differ diff --git a/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/runefieldsim_0_L0_0000.root b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/runefieldsim_0_L0_0000.root new file mode 100644 index 00000000..bb126c75 Binary files /dev/null and b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/runefieldsim_0_L0_0000.root differ diff --git a/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/runefieldsim_0_L1_0000.root b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/runefieldsim_0_L1_0000.root new file mode 100644 index 00000000..05583baa Binary files /dev/null and b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/runefieldsim_0_L1_0000.root differ diff --git a/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/runshowersim_0_L0_0000.root b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/runshowersim_0_L0_0000.root new file mode 100644 index 00000000..6edd8515 Binary files /dev/null and b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/runshowersim_0_L0_0000.root differ diff --git a/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/shower_5388-23832_L0_0000.root b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/shower_5388-23832_L0_0000.root new file mode 100644 index 00000000..4cfec111 Binary files /dev/null and b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/shower_5388-23832_L0_0000.root differ diff --git a/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/showersim_5388-23832_L0_0000.root b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/showersim_5388-23832_L0_0000.root new file mode 100644 index 00000000..dcbd0cbc Binary files /dev/null and b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/showersim_5388-23832_L0_0000.root differ diff --git a/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/voltage_5388-23832_L0_0000.root b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/voltage_5388-23832_L0_0000.root new file mode 100644 index 00000000..efc09438 Binary files /dev/null and b/sim2root/Common/sim_Xiaodushan_20221026_000000_RUN0_CD_ZHAireS_0000/voltage_5388-23832_L0_0000.root differ