diff --git a/all_info_fits.py b/all_info_fits.py index d5e6d884b95f7641fdb3a2c9e274e4fb586c6134..9967e679a6b8243a7091e751f1621b73f4e0cb27 100644 --- a/all_info_fits.py +++ b/all_info_fits.py @@ -5,7 +5,7 @@ import numpy as np np.set_printoptions(threshold=np.inf) # Percorso del file FITS -file_path = '/home/alfonso/Scrivania/THESEUS/xgis_m7-main/msk_55_4752313.fits' +file_path = '/home/alfonso/Scrivania/THESEUS/xgis_m7-main/fits_matrix/matrix.fits' # Apri il file FITS with fits.open(file_path) as hdul: diff --git a/primo.mac b/primo.mac index 10b70ff60710408cb8d03548037bfa97ad5c88a0..8843ae35406b692ee99f20871054a5603bf9fe6c 100644 --- a/primo.mac +++ b/primo.mac @@ -2,8 +2,8 @@ /gps/pos/type Plane /gps/pos/shape Square -/gps/pos/halfx 50 cm -/gps/pos/halfy 50 cm +/gps/pos/halfx 100 cm +/gps/pos/halfy 100 cm /gps/ang/type cos /process/em/fluo true @@ -15,5 +15,5 @@ /gps/pos/rot2 1.000 0.000 -0.000 /gps/direction -0.000 -0.000 -1.000 -/gps/energy 2.0 keV -/run/beamOn 100 +/gps/energy 154.0 keV +/run/beamOn 1000000 diff --git a/python/__pycache__/arf.cpython-312.pyc b/python/__pycache__/arf.cpython-312.pyc index 2db2af91a97d462c470049d34ba658ed5c399b49..2d3414f352cf78adcf19ad8419f13b93dbcf388b 100644 Binary files a/python/__pycache__/arf.cpython-312.pyc and b/python/__pycache__/arf.cpython-312.pyc differ diff --git a/python/__pycache__/convertRoot2Fits.cpython-38.pyc b/python/__pycache__/convertRoot2Fits.cpython-38.pyc index f6ed7417ae74281b26d13fd743b02b672a00a171..f5a944a5abe16fceea8d616c9f72a74aac424ce6 100644 Binary files a/python/__pycache__/convertRoot2Fits.cpython-38.pyc and b/python/__pycache__/convertRoot2Fits.cpython-38.pyc differ diff --git a/python/__pycache__/data.cpython-312.pyc b/python/__pycache__/data.cpython-312.pyc index e7c2441e7c75368b42ab6b4a3f46e7ae6ae9fef0..e564a2e66bb612919c0db9795d5e79ea546ac89a 100644 Binary files a/python/__pycache__/data.cpython-312.pyc and b/python/__pycache__/data.cpython-312.pyc differ diff --git a/python/__pycache__/macro.cpython-38.pyc b/python/__pycache__/macro.cpython-38.pyc index 565c67d7beb22aab33c68cfa61e710c397995b14..6e49a213eb95a5301658c941acc49f1b51da31c1 100644 Binary files a/python/__pycache__/macro.cpython-38.pyc and b/python/__pycache__/macro.cpython-38.pyc differ diff --git a/python/__pycache__/noisy.cpython-312.pyc b/python/__pycache__/noisy.cpython-312.pyc new file mode 100644 index 0000000000000000000000000000000000000000..ee34848c5d0ceccc919947ab44d62918b706ce9d Binary files /dev/null and b/python/__pycache__/noisy.cpython-312.pyc differ diff --git a/python/__pycache__/rmf.cpython-312.pyc b/python/__pycache__/rmf.cpython-312.pyc index 9598531da4f4748dd03cc387c8ad20056213d6ce..49914d38315c95fa98e449571bcdc9b7a5a8e487 100644 Binary files a/python/__pycache__/rmf.cpython-312.pyc and b/python/__pycache__/rmf.cpython-312.pyc differ diff --git a/python/aeFF.py b/python/aeFF.py index c5ea7bbdd01c926f8b947e3a00a08a48d82fe8cc..f8c2776584b150c79db179d59bb54f802fc261c8 100644 --- a/python/aeFF.py +++ b/python/aeFF.py @@ -12,7 +12,7 @@ def main(): # 1 Estract data from fits file and give these output file_dat file_dat = "areaEfficace.dat" data_estraction(file_dat) - + sys.exit() #2 Filter data data = np.loadtxt(file_dat, dtype=str) diff --git a/python/data.py b/python/data.py index 49c871b9485e68d4a060c0e3abf0064cc6b5de92..2c1b41ee970fd3090b710dad6b647d56492431c2 100644 --- a/python/data.py +++ b/python/data.py @@ -6,7 +6,7 @@ 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): +def write_fits_dep(energy, theta, phi, type, en_dep_noisy, en_dep, name_id, pos_id): # Write FITS file # "Null" primary array @@ -15,7 +15,8 @@ def write_fits_dep(energy, theta, phi, type, en_dep_noisy, en_dep): # Extension tbhdu = pyfits.BinTableHDU.from_columns([ pyfits.Column(name='En_dep_Noise', format='1D', unit='keV', array = en_dep_noisy), - pyfits.Column(name='En_dep', format='1D', unit='keV', array = en_dep) + pyfits.Column(name=name_id, format='1I', array = pos_id) + #pyfits.Column(name='En_dep', format='1D', unit='keV', array = en_dep) ]) tbhdu.header['EXTNAME'] = 'EVENTS' @@ -57,11 +58,14 @@ def get_data(files): with pyfits.open(files) as hdul: data = hdul[1].data tot_events = len(data) + pixel_id_col = data['Pixel_ID'] + pixel_id = pixel_id_col[pixel_id_col != 0] # ! dopo il run metti -1000 anche qui scint_id_col = data['Scint_ID'] + scint_id = scint_id_col[scint_id_col != -1000] energy = data['En_Primary'][0] en_dep = data['En_dep'] 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_S = tot_events - eventi_X ratio_X = eventi_X/tot_events @@ -70,16 +74,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 + 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 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 in results: - - X_fits = write_fits_dep(energy, theta, phi, 'X', X_en_dep_noisy, X_en_dep) - S_fits = write_fits_dep(energy, theta, phi, 'S', S_en_dep_noisy, S_en_dep) + 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: + + 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) f.write(f"{energy:.1f}\t{ratio_X:.3f}\t{ratio_S:.3f}\t{theta}\t{phi}\t{X_fits}\t{S_fits}\n") diff --git a/python/macro.py b/python/macro.py index 4c18daedd1c28c3d1917b4a54e72712d8636e9d8..3f46c52255bae83731781272d0e2375797fdb8ec 100644 --- a/python/macro.py +++ b/python/macro.py @@ -12,8 +12,8 @@ def create_macro(macro_name,energy,x,y,z,rot1_x,rot1_y,rot1_z,rot2_x,rot2_y,rot2 f.write("\n") f.write("/gps/pos/type Plane\n") f.write("/gps/pos/shape Square\n") - f.write("/gps/pos/halfx 50 cm\n") - f.write("/gps/pos/halfy 50 cm\n") + f.write("/gps/pos/halfx 100 cm\n") + f.write("/gps/pos/halfy 100 cm\n") f.write("/gps/ang/type cos\n") f.write("\n") f.write("/process/em/fluo true\n") @@ -26,7 +26,7 @@ def create_macro(macro_name,energy,x,y,z,rot1_x,rot1_y,rot1_z,rot2_x,rot2_y,rot2 f.write(f"/gps/direction {dir_x:.3f} {dir_y:.3f} {dir_z:.3f}\n") f.write("\n") f.write(f"/gps/energy {energy} keV\n") - f.write("/run/beamOn 10000\n") + f.write("/run/beamOn 1000000\n") print(f"File '{macro_name}' creato correttamente") except Exception as e: print(f"Errore durante la creazione del file: {e}") diff --git a/python/matrix.py b/python/matrix.py new file mode 100644 index 0000000000000000000000000000000000000000..ec1c547211ebf6108e96b0b91f589eeaec2126a6 --- /dev/null +++ b/python/matrix.py @@ -0,0 +1,72 @@ +import numpy as np +import astropy.io.fits as pyfits +import matplotlib.pyplot as plt +import os +np.set_printoptions(threshold=np.inf) + +file = "/home/alfonso/Scrivania/THESEUS/xgis_m7-main/fits_dep/X_en_dep_10.0_0.0_0.0.fits" +#file = "/home/alfonso/Scrivania/THESEUS/xgis_m7-main/fits_dep/S_en_dep_60.0_0.0_0.0.fits" + +with pyfits.open(file) as hdul: + data = hdul[1].data + pixel_id = data['Pixel_ID'] + #scint_id = data['Scint_ID'] + +events = pixel_id +#events = scint_id + +# matrice di zeri 89x89 +counts_matrix = np.zeros((89, 89), dtype=int) + +# per ogni eventoe aggiorna il conteggio nel pixel corrispondente +for event in events: + # dall'array trova le coo del pixel + row = event // 89 # divisione intera per riga + col = event % 89 # resto per colonna + counts_matrix[row, col] += 1 + +#print(counts_matrix) + +dead_matrix = np.ones((89, 89), dtype=int) + +# Imposta a zero ogni 9 righe e ogni 9 colonne +dead_matrix[8::9, :] = 0 # Ogni nona riga +dead_matrix[:, 8::9] = 0 # Ogni nona colonna + +IMMAGINE = counts_matrix*dead_matrix + +print(IMMAGINE) + +elements_sum = np.sum(IMMAGINE) +print("La somma totale dei valori nella matrice รจ:", elements_sum) + + +# Write FITS file +# "Null" primary array +prhdu = pyfits.PrimaryHDU() +imhdu = pyfits.ImageHDU(data=IMMAGINE) + + +imhdu.header['EXTNAME'] = 'MATRIX' +imhdu.header['TELESCOP'] = 'THESEUS' +imhdu.header['INSTRUME'] = 'XGIS' +imhdu.header['COMMENT'] = "Alfonso Pisapia" + +hdulist = pyfits.HDUList([prhdu, imhdu]) + +fits_matrix_dir = "/home/alfonso/Scrivania/THESEUS/xgis_m7-main/fits_matrix" +if not os.path.exists(fits_matrix_dir): + os.makedirs(fits_matrix_dir) + +file_name = "matrix.fits" + +output_file = os.path.join(fits_matrix_dir, file_name) +hdulist.writeto(output_file, overwrite=True) + + +plt.imshow(IMMAGINE, cmap='viridis', interpolation='nearest') +plt.colorbar(label='Conteggi') +plt.title('Matrice dei Conteggi degli Eventi') +plt.xlabel('Colonne') +plt.ylabel('Righe') +plt.show() \ No newline at end of file diff --git a/python/simulazioni.sh b/python/simulazioni.sh index f04090e74932cb93d28ecd84f53e38de914b370f..f0485b1b757de0f74c06a757d1e8599ad570dba0 100755 --- a/python/simulazioni.sh +++ b/python/simulazioni.sh @@ -2,11 +2,11 @@ # Iterazioni su energy, theta e phi, potrei chiedere in input lo step e Neventi for energy in $(seq 2 19 154); do # 2 2 152 - for theta in $(seq 0 30 90); do - for phi in $(seq 0 90 360); do - echo "Lancio simulazione con energy: ${energy} keV, Theta: ${theta}, Phi: ${phi}" - # Chiama lo script Python passando i parametri - python3.8 main.py $energy $theta $phi - done - done + # for theta in $(seq 0 30 90); do + # for phi in $(seq 0 90 360); do + # echo "Lancio simulazione con energy: ${energy} keV, Theta: ${theta}, Phi: ${phi}" + # # Chiama lo script Python passando i parametri + python3.8 main.py $energy 0 0 #$theta $phi + # done + # done done \ No newline at end of file diff --git a/src/UserRun.cc b/src/UserRun.cc index 5e5f5332d1116e00cc5a00d3b21890262d44867b..4b463c5a6cbf4b457b25b0b01fdf42076128c6fa 100644 --- a/src/UserRun.cc +++ b/src/UserRun.cc @@ -138,7 +138,7 @@ void UserRun::RecordEvent(const G4Event* event) analysisManager->FillNtupleDColumn(0, eventID); analysisManager->FillNtupleDColumn(1, energyDeposit/keV); analysisManager->FillNtupleDColumn(2, xpixel); // ScintID - analysisManager->FillNtupleDColumn(3, ypixel); // PixelID is the slice number for scintillator + analysisManager->FillNtupleDColumn(3, -1000); // ypixel -> PixelID is the slice number for scintillator //! per ora -1000 PixelID is undefined for Scintillator analysisManager->FillNtupleDColumn(4, locpixelX/cm); analysisManager->FillNtupleDColumn(5, locpixelY/cm); analysisManager->FillNtupleDColumn(6, locpixelZ/cm);