Skip to content
Snippets Groups Projects
Commit 58fc7424 authored by Alfonso's avatar Alfonso
Browse files

per ritrovare kapton

parent 0002b01d
Branches
Tags
No related merge requests found
...@@ -12,7 +12,7 @@ from rmf import rmf_file ...@@ -12,7 +12,7 @@ from rmf import rmf_file
def main(): def main():
# 1 Estract data from fits file and give this output file_dat # 1 Estract data from fits file and give this output file_dat
file_dat = "areaEfficace.dat" file_dat = "areaEfficace.dat" # è solo un file di testo che mi lista i dati, serve per scrivere rmf e arf, passarci sempre anche se non faccio rmf e arf
data_estraction(file_dat) data_estraction(file_dat)
#2 Filter data #2 Filter data
......
...@@ -6,7 +6,9 @@ import astropy.io.fits as pyfits ...@@ -6,7 +6,9 @@ import astropy.io.fits as pyfits
from noisy import add_noise_X, add_noise_S from noisy import add_noise_X, add_noise_S
def write_fits_dep(energy, theta, phi, type, en_dep_noisy, en_dep, name_id, pos_id): def write_fits_dep(energy, en_primary, theta, phi, type, en_dep_noisy, en_dep, name_id, pos_id, x_pixel, y_pixel): # ! è qui che decido quali informazioni sulla simulazione voglio nei file fits_dep
# Write FITS file # Write FITS file
# "Null" primary array # "Null" primary array
...@@ -14,18 +16,19 @@ def write_fits_dep(energy, theta, phi, type, en_dep_noisy, en_dep, name_id, pos_ ...@@ -14,18 +16,19 @@ def write_fits_dep(energy, theta, phi, type, en_dep_noisy, en_dep, name_id, pos_
# Extension # Extension
tbhdu = pyfits.BinTableHDU.from_columns([ tbhdu = pyfits.BinTableHDU.from_columns([
pyfits.Column(name='En_dep_Noise', format='1D', unit='keV', array = en_dep_noisy), #pyfits.Column(name='En_dep_Noise', format='1D', unit='keV', array = en_dep_noisy),
pyfits.Column(name=name_id, format='1I', array = pos_id) pyfits.Column(name='En_Primary', format='1D', unit='keV', array = en_primary),
#pyfits.Column(name='(X,Y)', format='2I', array = (x_pixel,y_pixel))
#pyfits.Column(name='En_dep', format='1D', unit='keV', array = en_dep) #pyfits.Column(name='En_dep', format='1D', unit='keV', array = en_dep)
]) ])
tbhdu.header['EXTNAME'] = 'EVENTS' tbhdu.header['EXTNAME'] = 'EVENTS'
tbhdu.header['TELESCOP'] = 'THESEUS' tbhdu.header['TELESCOP'] = 'THESEUS'
tbhdu.header['INSTRUME'] = 'XGIS' tbhdu.header['INSTRUME'] = 'XGIS'
tbhdu.header['COMMENT'] = "Simulated energy deposits with noise" tbhdu.header['COMMENT'] = "Simulated GRB spectrum, Band funcition alpha1 = -1, alpha2 = -3, Ebreak = 100 "
tbhdu.header['TUNIT1'] = "keV" tbhdu.header['TUNIT1'] = "keV"
tbhdu.header['TYPE'] = type tbhdu.header['TYPE'] = type
tbhdu.header['ENERGY'] = energy #tbhdu.header['ENERGY'] = energy
tbhdu.header['THETA'] = theta tbhdu.header['THETA'] = theta
tbhdu.header['PHI'] = phi tbhdu.header['PHI'] = phi
tbhdu.header['COMMENT'] = "Alfonso Pisapia" tbhdu.header['COMMENT'] = "Alfonso Pisapia"
...@@ -36,14 +39,13 @@ def write_fits_dep(energy, theta, phi, type, en_dep_noisy, en_dep, name_id, pos_ ...@@ -36,14 +39,13 @@ def write_fits_dep(energy, theta, phi, type, en_dep_noisy, en_dep, name_id, pos_
if not os.path.exists(fits_dep_dir): if not os.path.exists(fits_dep_dir):
os.makedirs(fits_dep_dir) os.makedirs(fits_dep_dir)
file_name = f"{type}_en_dep_{energy:.1f}_{theta}_{phi}.fits" file_name = f"{type}_grb_{theta}_{phi}.fits"
output_file = os.path.join(fits_dep_dir, file_name) output_file = os.path.join(fits_dep_dir, file_name)
hdulist.writeto(output_file, overwrite=True) hdulist.writeto(output_file, overwrite=True)
return file_name return file_name
# Funzione che prende in input un singolo file FITS da processare. # Funzione che prende in input un singolo file FITS da processare.
def get_data(files): def get_data(files):
# Estrazione di energia, theta e phi dal nome del file # Estrazione di energia, theta e phi dal nome del file
...@@ -54,7 +56,7 @@ def get_data(files): ...@@ -54,7 +56,7 @@ def get_data(files):
theta = float(name_parts[2]) theta = float(name_parts[2])
phi = float(name_parts[3]) phi = float(name_parts[3])
# Apertura del file fits e lettura dei dati # Apertura del file fits e lettura dei dati # ! in base alle info che voglio nel fits_dep qui devo richiamare le info
with pyfits.open(files) as hdul: with pyfits.open(files) as hdul:
data = hdul[1].data data = hdul[1].data
tot_events = len(data) tot_events = len(data)
...@@ -63,7 +65,10 @@ def get_data(files): ...@@ -63,7 +65,10 @@ def get_data(files):
scint_id_col = data['Scint_ID'] scint_id_col = data['Scint_ID']
scint_id = scint_id_col[scint_id_col != -1000] scint_id = scint_id_col[scint_id_col != -1000]
energy = data['En_Primary'][0] energy = data['En_Primary'][0]
en_primary = data['En_Primary']
en_dep = data['En_dep'] en_dep = data['En_dep']
x_pixel = data['X_Pixel']
y_pixel = data['Y_Pixel']
X_en_dep = en_dep[scint_id_col == -1000] X_en_dep = en_dep[scint_id_col == -1000]
S_en_dep = en_dep[scint_id_col != -1000] S_en_dep = en_dep[scint_id_col != -1000]
eventi_X = (scint_id_col == -1000).sum() eventi_X = (scint_id_col == -1000).sum()
...@@ -74,16 +79,16 @@ def get_data(files): ...@@ -74,16 +79,16 @@ def get_data(files):
X_en_dep_noisy = [add_noise_X(X_dep) for X_dep in X_en_dep] # list comprehension di Python, una struttura compatta che consente di applicare una funzione a ciascun elemento di un array o di una lista, creando direttamente un nuovo array con i risultati. X_en_dep_noisy = [add_noise_X(X_dep) for X_dep in X_en_dep] # list comprehension di Python, una struttura compatta che consente di applicare una funzione a ciascun elemento di un array o di una lista, creando direttamente un nuovo array con i risultati.
S_en_dep_noisy = [add_noise_S(S_dep) for S_dep in S_en_dep] S_en_dep_noisy = [add_noise_S(S_dep) for S_dep in S_en_dep]
return energy, ratio_X, ratio_S, theta, phi, X_en_dep_noisy, S_en_dep_noisy, X_en_dep, S_en_dep, pixel_id, scint_id return energy, en_primary, ratio_X, ratio_S, theta, phi, X_en_dep_noisy, S_en_dep_noisy, X_en_dep, S_en_dep, pixel_id, scint_id, x_pixel, y_pixel
# scrive sia il file di testo che i file fits
def write_results(file_dat, results): def write_results(file_dat, results):
with open(file_dat, "w") as f: with open(file_dat, "w") as f:
for energy, ratio_X, ratio_S, theta, phi, X_en_dep_noisy, S_en_dep_noisy, X_en_dep, S_en_dep, pixel_id, scint_id in results: for energy, en_primary, ratio_X, ratio_S, theta, phi, X_en_dep_noisy, S_en_dep_noisy, X_en_dep, S_en_dep, pixel_id, scint_id in results:
X_fits = write_fits_dep(energy, theta, phi, 'X', X_en_dep_noisy, X_en_dep, 'Pixel_ID', pixel_id) X_fits = write_fits_dep(energy, en_primary, theta, phi, 'X', X_en_dep_noisy, X_en_dep, 'Pixel_ID', pixel_id)
S_fits = write_fits_dep(energy, theta, phi, 'S', S_en_dep_noisy, S_en_dep, 'Scint_ID', scint_id) S_fits = write_fits_dep(energy, en_primary, theta, phi, 'S', S_en_dep_noisy, S_en_dep, 'Scint_ID', scint_id)
f.write(f"{energy:.1f}\t{ratio_X:.3f}\t{ratio_S:.3f}\t{theta}\t{phi}\t{X_fits}\t{S_fits}\n") f.write(f"{energy:.1f}\t{ratio_X:.3f}\t{ratio_S:.3f}\t{theta}\t{phi}\t{X_fits}\t{S_fits}\n")
...@@ -103,4 +108,4 @@ def data_estraction(file_dat): ...@@ -103,4 +108,4 @@ def data_estraction(file_dat):
with Pool() as pool: with Pool() as pool:
results = pool.map(get_data, files) # qui files è l'argomento che passi a get_data results = pool.map(get_data, files) # qui files è l'argomento che passi a get_data
write_results(file_dat, results) write_results(file_dat, results) # per tutti i file fits che trova, in parallelo, ottiene i dati e li scrive in file_dat
\ No newline at end of file \ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment