diff --git a/python/aeFF.py b/python/aeFF.py index 01b8aa4e5ae2082e844f6efb430b42e43fcc3a95..aa858852415e91a9a4d961cfce578cdbba48e66b 100644 --- a/python/aeFF.py +++ b/python/aeFF.py @@ -12,7 +12,7 @@ from rmf import rmf_file def main(): # 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) #2 Filter data diff --git a/python/convertRoot2Fits.py b/python/convertRoot2Fits.py index abe02c5b51fbee0120d31b7bad6dee5403f7af33..090dbff732a1e834d771e5f546182b30de394f29 100644 --- a/python/convertRoot2Fits.py +++ b/python/convertRoot2Fits.py @@ -65,7 +65,7 @@ def generate_fits(energy, theta_deg, phi_deg): pyfits.Column(name='Pixel_ID', format='1K', array=t_SL[sorting]), pyfits.Column(name='X_Pixel', format='1D', unit='cm', array=t_XH[sorting]), - + pyfits.Column(name='Y_Pixel', format='1D', unit='cm', array=t_YH[sorting]), pyfits.Column(name='Z_Pixel', format='1D', unit='cm', array=t_ZH[sorting]), diff --git a/python/data.py b/python/data.py index a403951b9e25ed4006514434ec7e95631658ebd4..af7f0a61d383207ae9887290e277a0fa2fcb6cab 100644 --- a/python/data.py +++ b/python/data.py @@ -6,7 +6,9 @@ import astropy.io.fits as pyfits 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 # "Null" primary array @@ -14,18 +16,19 @@ def write_fits_dep(energy, theta, phi, type, en_dep_noisy, en_dep, name_id, pos_ # Extension tbhdu = pyfits.BinTableHDU.from_columns([ - 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_dep_Noise', format='1D', unit='keV', array = en_dep_noisy), + 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) ]) tbhdu.header['EXTNAME'] = 'EVENTS' tbhdu.header['TELESCOP'] = 'THESEUS' 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['TYPE'] = type - tbhdu.header['ENERGY'] = energy + #tbhdu.header['ENERGY'] = energy tbhdu.header['THETA'] = theta tbhdu.header['PHI'] = phi tbhdu.header['COMMENT'] = "Alfonso Pisapia" @@ -36,13 +39,12 @@ def write_fits_dep(energy, theta, phi, type, en_dep_noisy, en_dep, name_id, pos_ if not os.path.exists(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) hdulist.writeto(output_file, overwrite=True) return file_name - # Funzione che prende in input un singolo file FITS da processare. def get_data(files): @@ -54,7 +56,7 @@ def get_data(files): theta = float(name_parts[2]) 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: data = hdul[1].data tot_events = len(data) @@ -63,7 +65,10 @@ def get_data(files): scint_id_col = data['Scint_ID'] scint_id = scint_id_col[scint_id_col != -1000] energy = data['En_Primary'][0] + en_primary = data['En_Primary'] en_dep = data['En_dep'] + x_pixel = data['X_Pixel'] + y_pixel = data['Y_Pixel'] X_en_dep = en_dep[scint_id_col == -1000] S_en_dep = en_dep[scint_id_col != -1000] eventi_X = (scint_id_col == -1000).sum() @@ -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. 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): 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) - S_fits = write_fits_dep(energy, theta, phi, 'S', S_en_dep_noisy, S_en_dep, 'Scint_ID', scint_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, 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") @@ -103,4 +108,4 @@ def data_estraction(file_dat): with Pool() as pool: results = pool.map(get_data, files) # qui files è l'argomento che passi a get_data - write_results(file_dat, results) \ No newline at end of file + 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