diff --git a/.vscode/settings.json b/.vscode/settings.json
new file mode 100644
index 0000000000000000000000000000000000000000..7c3b1b378ab8d7274450b91ae6e247ddf725121f
--- /dev/null
+++ b/.vscode/settings.json
@@ -0,0 +1,49 @@
+{
+ "files.associations": {
+ "array": "cpp",
+ "atomic": "cpp",
+ "bit": "cpp",
+ "*.tcc": "cpp",
+ "cctype": "cpp",
+ "clocale": "cpp",
+ "cmath": "cpp",
+ "cstdarg": "cpp",
+ "cstddef": "cpp",
+ "cstdint": "cpp",
+ "cstdio": "cpp",
+ "cstdlib": "cpp",
+ "cwchar": "cpp",
+ "cwctype": "cpp",
+ "deque": "cpp",
+ "map": "cpp",
+ "unordered_map": "cpp",
+ "vector": "cpp",
+ "exception": "cpp",
+ "algorithm": "cpp",
+ "functional": "cpp",
+ "iterator": "cpp",
+ "memory": "cpp",
+ "memory_resource": "cpp",
+ "numeric": "cpp",
+ "optional": "cpp",
+ "random": "cpp",
+ "string": "cpp",
+ "string_view": "cpp",
+ "system_error": "cpp",
+ "tuple": "cpp",
+ "type_traits": "cpp",
+ "utility": "cpp",
+ "fstream": "cpp",
+ "initializer_list": "cpp",
+ "iosfwd": "cpp",
+ "iostream": "cpp",
+ "istream": "cpp",
+ "limits": "cpp",
+ "new": "cpp",
+ "ostream": "cpp",
+ "sstream": "cpp",
+ "stdexcept": "cpp",
+ "streambuf": "cpp",
+ "typeinfo": "cpp"
+ }
+}
\ No newline at end of file
diff --git a/all_info_fits.py b/all_info_fits.py
index a1a9e7a093e0ee3b0017afab8a9c588aa00d221c..d5e6d884b95f7641fdb3a2c9e274e4fb586c6134 100644
--- a/all_info_fits.py
+++ b/all_info_fits.py
@@ -1,10 +1,11 @@
-#printa tutte le info di un file fits
-#sostituisci il path del file .fits
-
from astropy.io import fits
+import numpy as np
+
+# Configura numpy per visualizzare tutti i dati senza troncamenti
+np.set_printoptions(threshold=np.inf)
# Percorso del file FITS
-file_path = '/home/alfonso/Scrivania/THESEUS/XGIS-X_0deg_reqEOL_v8.rmf'
+file_path = '/home/alfonso/Scrivania/THESEUS/xgis_m7-main/msk_55_4752313.fits'
# Apri il file FITS
with fits.open(file_path) as hdul:
@@ -12,8 +13,9 @@ with fits.open(file_path) as hdul:
hdul.info()
# Itera su tutte le HDU e stampa l'header e i dati
- for hdu in hdul:
- print(f"\nHeader HDU {hdul.index(hdu)}:")
+ for idx, hdu in enumerate(hdul):
+ print(f"\nHeader HDU {idx}:")
print(hdu.header)
print("\nDati:")
- print(hdu.data)
+ print(hdu.data) # Stampa i dati completi senza troncamenti
+
diff --git a/detector_params.inp b/detector_params.inp
index 87b086ad2686f0086ce4264fc3861e4896611dba..9f564f66bfcf518ebca722a19fc75048dc1cab0d 100644
--- a/detector_params.inp
+++ b/detector_params.inp
@@ -11,9 +11,9 @@ BUS_THICK = 20. # Bottom bus dummy thickness (in cm)
BUS_HEIGHT = 150. # Spacecraft bus height (in cm)
BUS_SIDE = 100. # Spacecraft bus side (in cm)
SDD_THICK = 450.0 # SDD thickness (in um)
-MASK_X = 56.1 # X side of the coded mask (in cm)
-MASK_Y = 56.1 # Y side of the coded mask (in cm)
-MASK_THICK = 1000. # Thickness of the coded mask (in um)
+MASK_X = 56.5 # X side of the coded mask (in cm) # 56.1 old
+MASK_Y = 56.5 # Y side of the coded mask (in cm)
+MASK_THICK = 500. # Thickness of the coded mask (in um) # 1000 old
MASK_FOCAL = 63. # Focal distance btw. mask and detector (in cm)
COLL_UPPER_1_THICK = 1000.0 # Upper internal collimator thickness (in um)
COLL_UPPER_2_THICK = 250.0 # Upper external collimator thickness (in um)
@@ -23,7 +23,7 @@ OPTICAL_COUPLER_THICK = 1.0 # Optical coupler thickness (in mm)
# Module parameters
N_SIDE = 89 # Number of scintillator bars per side (89x89=7921-81)
N_SLICES_PER_BAR = 1 # Number of voxels (slices) for each scintillator bar
-N_MASK_ELEM_PER_SIDE = 16 # Number of mask elements per side (+16 vuoti)
+N_MASK_ELEM_PER_SIDE = 55 # Number of mask elements per side #16 old
# Materials
diff --git a/include/MaskParameterisation.hh b/include/MaskParameterisation.hh
index 0816daee6d4387528be57ab62bdcbb390d990df8..46b5df6256d0f969c280d9a6b36a0f7a53270d69 100644
--- a/include/MaskParameterisation.hh
+++ b/include/MaskParameterisation.hh
@@ -4,6 +4,11 @@
#include "globals.hh"
#include "G4VPVParameterisation.hh"
+#include "G4ThreeVector.hh"
+#include <vector>
+#include <utility>
+#include <string>
+
class G4VPhysicalVolume;
class G4Box;
@@ -25,12 +30,13 @@ class G4BREPSolidPolyhedra;
class MaskParameterisation : public G4VPVParameterisation
{
public:
- MaskParameterisation(G4double, G4double, G4int);
+ MaskParameterisation(G4double side_length_X, G4double side_length_Y, G4int n_side, G4double gap, const std::string& filename);
virtual ~MaskParameterisation();
- void ComputeTransformation (const G4int copyNo, G4VPhysicalVolume* physVol) const;
-
+ void ComputeTransformation (const G4int copyNo, G4VPhysicalVolume* physVol) const override;
+
+ int GetTotalZeros() const;
private: // Dummy declarations to get rid of warnings...
void ComputeDimensions (G4Trd&,const G4int,const G4VPhysicalVolume*) const {}
@@ -50,7 +56,12 @@ class MaskParameterisation : public G4VPVParameterisation
private:
G4double fSideX;
G4double fSideY;
+ G4double fGap;
G4int nSide;
+ std::vector<std::pair<int, int>> zeroPositions;
+
+ int LoadZeroPositions(const std::string& filename);
+
};
diff --git a/python/a.out b/python/a.out
new file mode 100755
index 0000000000000000000000000000000000000000..67b98a4d229102453ed94ee2a9591fbce78cb5a2
Binary files /dev/null and b/python/a.out differ
diff --git a/python/aeFF.py b/python/aeFF.py
index e07523ee65be4db676a36ceb7a56dcc690010b7e..2f800aeaa5ae76f1b003477770f586146d214bb2 100644
--- a/python/aeFF.py
+++ b/python/aeFF.py
@@ -1,4 +1,5 @@
# come faccio a sapere se il mio codice è robusto
+# aEFF deve essere corretta nei dati
import numpy as np
import sys
@@ -44,18 +45,15 @@ def main():
S_en_dep_file_filtered = S_en_dep_file[mask]
# Salva i dati per il plot
- plot_data_X.append((energy_filtered, ratio_X_filtered, theta_val, phi_val))
- plot_data_S.append((energy_filtered, ratio_S_filtered, theta_val, phi_val))
+ #plot_data_X.append((energy_filtered, ratio_X_filtered, theta_val, phi_val))
+ #plot_data_S.append((energy_filtered, ratio_S_filtered, theta_val, phi_val))
# 3 ARF
#arf_file(energy_filtered, ratio_X_filtered, ratio_S_filtered, theta_val, phi_val)
# 4 RMF
- #print(energy_filtered)
- #print(theta_val)
- #print(phi_val)
- rmf_file(energy_filtered, 'X', X_en_dep_file_filtered, theta_val, phi_val)
- rmf_file(energy_filtered, 'S', S_en_dep_file_filtered, theta_val, phi_val)
+ #rmf_file(energy_filtered, 'X', X_en_dep_file_filtered, theta_val, phi_val)
+ #rmf_file(energy_filtered, 'S', S_en_dep_file_filtered, theta_val, phi_val)
# approcci diversi per arf e rmf
# in arf essendo più semplice la chiamo una volta e dentro scrivo due volte le cose che deve fare
diff --git a/python/areaEfficace.arf b/python/areaEfficace.arf
deleted file mode 100644
index 16f5c4044f76926fd5a7af88716f011a1521cd18..0000000000000000000000000000000000000000
Binary files a/python/areaEfficace.arf and /dev/null differ
diff --git a/python/convertFits2Ascii.py b/python/convertFits2Ascii.py
new file mode 100644
index 0000000000000000000000000000000000000000..0d33c62affe30998ecde6bab9c598c5204456005
--- /dev/null
+++ b/python/convertFits2Ascii.py
@@ -0,0 +1,30 @@
+import astropy.io.fits as fits
+import numpy as np
+
+# ! posso anche lanciare questo script dal bash iniziale e gli do come argomento il nome del file
+
+mask_file_fits = "msk_55_4752313.fits"
+
+def fits_to_ascii(mask_file_fits, mask_file_ascii):
+
+ with fits.open(mask_file_fits) as hdul:
+ data = hdul[0].data
+
+ # Estrarre l'HDU primaria (dati e header)
+ primary_hdu = hdul[0]
+ data = primary_hdu.data
+ header = hdul[0].header
+
+ # Scrivi l'header e i dati in un file ASCII
+ with open(mask_file_ascii, 'w') as f:
+ # Scrivi l'header come commento
+ f.write("# Header Information\n")
+ for key, value in header.items():
+ f.write(f"# {key} = {value}\n")
+ f.write("# Data 0 = pieno 1 = vuoto\n")
+
+ # Scrivi i dati matrice nel file ASCII
+ np.savetxt(f, data, fmt="%d", delimiter=" ")
+
+# Esegui la funzione specificando il file FITS di input e il file ASCII di output
+fits_to_ascii('/home/alfonso/Scrivania/THESEUS/xgis_m7-main/python/msk_55_4752313.fits', 'mask.dat')
\ No newline at end of file
diff --git a/python/data.py b/python/data.py
index bf3e17e4ca2f0e7374402391062e9e1a652b19de..3f0c28217a745c24fc8bf2ba209b0ef05b371710 100644
--- a/python/data.py
+++ b/python/data.py
@@ -3,6 +3,9 @@ import numpy as np
from multiprocessing import Pool
import astropy.io.fits as fits
+from noisy import add_noise_X, add_noise_S
+
+
# Funzione che prende in input un singolo file FITS da processare.
def get_data(files):
# Estrazione di energia, theta e phi dal nome del file
@@ -27,12 +30,19 @@ def get_data(files):
ratio_X = eventi_X/tot_events
ratio_S = eventi_S/tot_events
+ X_en_dep_noisy = [add_noise_X(X_dep) for X_dep in X_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, S_en_dep
def write_results(file_dat, results):
+
npy_dir = "/home/alfonso/Scrivania/THESEUS/xgis_m7-main/python/npy"
+ npy_noise_dir = "/home/alfonso/Scrivania/THESEUS/xgis_m7-main/python/npy_noise"
if not os.path.exists(npy_dir):
os.makedirs(npy_dir)
+ if not os.path.exists(npy_noise_dir):
+ os.makedirs(npy_noise_dir)
with open(file_dat, "w") as f:
for energy, ratio_X, ratio_S, theta, phi, X_en_dep, S_en_dep in results:
@@ -48,6 +58,14 @@ def write_results(file_dat, results):
f.write(f"{energy:.1f}\t{ratio_X:.3f}\t{ratio_S:.3f}\t{theta}\t{phi}\t{X_npy_file}\t{S_npy_file}\n")
+ # Add noise to Energy Depoists Files and save in another folder
+ npy_noise_dir = "/home/alfonso/Scrivania/THESEUS/xgis_m7-main/python/npy_noise"
+ if not os.path.exists(npy_noise_dir):
+ os.makedirs(npy_noise_dir)
+
+ add_noise_X()
+ add_noise_S()
+
# Estrae i dati, processa in parallelo tutti i file .fits presenti nella folder, usa le due funzioni già definite
def data_estraction(file_dat):
# Crea una lista di file FITS trovati nella folder. Sorted per leggerli in ordine
diff --git a/python/noisy.py b/python/noisy.py
new file mode 100644
index 0000000000000000000000000000000000000000..7e70bfe1f8f96ff870ff52b4dc4f93058445380c
--- /dev/null
+++ b/python/noisy.py
@@ -0,0 +1,85 @@
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def add_noise_X(input):
+
+
+
+ w = 3.66 # eV/e-
+ fano_factor = 0.125 # for semiconductor detectors
+ sigma_fano = np.sqrt(fano_factor*input/w) # keV*(e-)/eV -> e-rms
+ sigma_el = 30.0 # e- rms in termini di devizaioni standard
+
+ variance = sigma_fano**2 + sigma_el**2 # Varianza della gaussiana
+ std_row = np.sqrt(variance)
+ std = (std_row*3.66)/1000 # (e-rms)*eV/(e-rms) --> keV
+
+ # Estrai un valore campionato dalla distribuzione gaussiana
+ output = np.random.normal(loc=input, scale=std)
+ #print(std)
+ #print("Valore estratto dalla distribuzione gaussiana:", output)
+
+
+ os.chdir("/home/alfonso/Scrivania/THESEUS/xgis_m7-main/python/npy")
+ array = np.load("X_en_dep_2.0_0.0_0.0.npy")
+ #print(array)
+
+ dep = [add_noise(input) for input in array] # 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_npy_noise_file = "X_en_dep_noy_2.0_0.0_0.0.npy"
+ X_en_dep_noise_file = os.path.join(npy_noise_dir, X_npy_noise_file)
+
+ new = np.save(X_en_dep_noise_file, dep)
+
+ loaded_dep = np.load(X_en_dep_noise_file)
+ print("Contenuto del nuovo file .npy con rumore aggiunto:")
+ print(loaded_dep)
+
+
+
+
+
+
+def add_noise_S(input):
+
+ light_efficiency = 25 # e-/keV
+ sigma_stat = input*light_efficiency
+ sigma_el = 30.0 # conto due volte perche ho due sdd
+
+ variance = sigma_stat**2 + 2*sigma_el**2
+
+ std = np.sqrt(variance)
+
+ # Estrai un valore campionato dalla distribuzione gaussiana
+ output = np.random.normal(loc=input, scale=std)
+ #print(std)
+ #print("Valore estratto dalla distribuzione gaussiana:", output)
+
+
+ return output
+
+
+
+
+
+
+ # # PLOT
+ # # Genera i valori per la curva gaussiana
+ # x = np.linspace(input - 4*std, input + 4*std, 1000)
+ # y = (1 / (std * np.sqrt(2 * np.pi))) * np.exp(-0.5 * ((x - input) / std)**2)
+ # # Plot della distribuzione gaussiana
+ # plt.plot(x, y, label="Distribuzione Gaussiana", color="blue")
+ # plt.axvline(output, color="red", linestyle="--", label=f"Valore campionato: {output:.2f}")
+ # # Aggiungi etichette e titolo
+ # plt.xlabel("Valori")
+ # plt.ylabel("Densità di probabilità")
+ # plt.title("Distribuzione Gaussiana con valore campionato evidenziato")
+ # plt.legend()
+ # plt.grid()
+ # # Mostra il grafico
+ # plt.show()
diff --git a/python/rmf.py b/python/rmf.py
index c064e45c613c6e23a3048b05aa11bb787b0a3deb..810bee27c03d7dda71e8e432fdc5647436f1460e 100644
--- a/python/rmf.py
+++ b/python/rmf.py
@@ -3,50 +3,49 @@ import os
import sys
import numpy as np
-
-def create_ebounds_header(extension):
- extension.header['EXTNAME'] = 'EBOUNDS'
+
+def create_matrix_header(extension):
+ extension.header['EXTNAME'] = 'MATRIX'
extension.header['TELESCOP'] = 'THESEUS'
extension.header['INSTRUME'] = 'XGIS'
extension.header['FILTER'] = 'NONE'
extension.header['CHANTYPE'] = 'PHA' # PI
- extension.header['DETCHANS'] = 'NONE' # the total number of raw detector PHA channels in the full (uncompressed) matrix
+ extension.header['DETCHANS'] = '31' # = bins, # the total number of raw detector PHA channels in the full (uncompressed) matrix.
extension.header['HDUCLASS'] = 'OGIP'
extension.header['HDUCLAS1'] = 'RESPONSE'
- extension.header['HDUCLAS2'] = 'EBOUNDS'
- extension.header['HDUVERS'] = '1.2.0'
+ extension.header['HDUCLAS2'] = 'RSP_MATRIX'
+ extension.header['HDUVERS'] = '1.3.0'
+ #extension.header['TLMIN#'] = 'NONE'# the first channel in the response. # is the column number for the F_CHAN column (see below).
# optional
+ #extension.header['NUMGRP'] = 'NONE' # the total number of channel subsets. The sum of the N_GRP column.
+ #extension.header['NUMELT'] = 'NONE' # the total number of response elements. The sum of the N_CHAN column
#extension.header['PHAFILE'] = 'NONE' # name of PHA file for which this file was produced
+ extension.header['LO_THRES '] = '0.0' # minimum probability threshold used to construct the matrix (matrix elements below this value are considered to zero and are not stored)
+ extension.header['HDUCLAS3'] = 'REDIST' # giving further details of the stored matrix Allowed values are: 'REDIST' for a matrix whose elements represent probabilities associated with the photon redistribution process only 'DETECTOR' for a matrix whose elements have been multiplied by all energy-dependent effects associated with detector (eg detector efficiency, window transmission etc). 'FULL' for a matrix whose elements have been multiplied by all energy-dependent effects associated with detector, optics, collimator, filters etc.
# obsolete (for old software)
extension.header['RMFVERSN'] = '1992A'
extension.header['HDUVERS1'] = '1.0.0'
- extension.header['HDUVERS2'] = '1.1.0'
-
-
-def create_matrix_header(extension):
- extension.header['EXTNAME'] = 'SPECRESP MATRIX'
+ extension.header['HDUVERS2'] = '1.3.0'
+
+
+def create_ebounds_header(extension):
+ extension.header['EXTNAME'] = 'EBOUNDS'
extension.header['TELESCOP'] = 'THESEUS'
extension.header['INSTRUME'] = 'XGIS'
extension.header['FILTER'] = 'NONE'
extension.header['CHANTYPE'] = 'PHA' # PI
- extension.header['DETCHANS'] = 'NONE' # the total number of raw detector PHA channels in the full (uncompressed) matrix.
+ extension.header['DETCHANS'] = '31' # = bins, the total number of raw detector PHA channels in the full (uncompressed) matrix
extension.header['HDUCLASS'] = 'OGIP'
extension.header['HDUCLAS1'] = 'RESPONSE'
- extension.header['HDUCLAS2'] = 'RSP_MATRIX'
- extension.header['HDUVERS'] = '1.3.0'
- #extension.header['TLMIN#'] = 'NONE'# the first channel in the response. # is the column number for the F_CHAN column (see below).
+ extension.header['HDUCLAS2'] = 'EBOUNDS'
+ extension.header['HDUVERS'] = '1.2.0'
# optional
- #extension.header['NUMGRP'] = 'NONE' # the total number of channel subsets. The sum of the N_GRP column.
- #extension.header['NUMELT'] = 'NONE' # the total number of response elements. The sum of the N_CHAN column
#extension.header['PHAFILE'] = 'NONE' # name of PHA file for which this file was produced
- extension.header['LO_THRES '] = '0' # minimum probability threshold used to construct the matrix (matrix elements below this value are considered to zero and are not stored)
- extension.header['HDUCLAS3'] = 'REDIST' # giving further details of the stored matrix Allowed values are: 'REDIST' for a matrix whose elements represent probabilities associated with the photon redistribution process only 'DETECTOR' for a matrix whose elements have been multiplied by all energy-dependent effects associated with detector (eg detector efficiency, window transmission etc). 'FULL' for a matrix whose elements have been multiplied by all energy-dependent effects associated with detector, optics, collimator, filters etc.
# obsolete (for old software)
extension.header['RMFVERSN'] = '1992A'
extension.header['HDUVERS1'] = '1.0.0'
- extension.header['HDUVERS2'] = '1.3.0'
-
- # altre mandatory che non ho capito
+ extension.header['HDUVERS2'] = '1.1.0'
+
def estract_row(en_dep_file, bins):
@@ -87,6 +86,7 @@ def estract_row(en_dep_file, bins):
return n_grp, f_chan, n_chan, counts_norm # row_matrix = counts_norm
+
def rmf_file(energy, type, en_dep_file_filt, theta, phi):
try:
os.chdir("/home/alfonso/Scrivania/THESEUS/xgis_m7-main/python/npy")
@@ -121,14 +121,6 @@ def rmf_file(energy, type, en_dep_file_filt, theta, phi):
# Creazione di "Null" primary array
primary_hdu = pyfits.PrimaryHDU()
- # Creazione di ebounds BinTableHDU
- ebounds_bin_table_hdu = pyfits.BinTableHDU.from_columns([
- pyfits.Column(name='CHANNEL', format='1I', array = channels),
- pyfits.Column(name='E_MIN', format='1E', unit='keV', array = bins[:-1]),
- pyfits.Column(name='E_MAX', format='1E', unit='keV', array = bins[1:]),
- ])
- create_ebounds_header(ebounds_bin_table_hdu)
-
# Creazione di matrice BinTableHDU
matrix_bin_table_hdu = pyfits.BinTableHDU.from_columns([
pyfits.Column(name='ENERG_LO', format='1E', unit='keV', array = energy[:-1]),
@@ -140,8 +132,16 @@ def rmf_file(energy, type, en_dep_file_filt, theta, phi):
])
create_matrix_header(matrix_bin_table_hdu)
+ # Creazione di ebounds BinTableHDU
+ ebounds_bin_table_hdu = pyfits.BinTableHDU.from_columns([
+ pyfits.Column(name='CHANNEL', format='1I', array = channels),
+ pyfits.Column(name='E_MIN', format='1E', unit='keV', array = bins[:-1]),
+ pyfits.Column(name='E_MAX', format='1E', unit='keV', array = bins[1:]),
+ ])
+ create_ebounds_header(ebounds_bin_table_hdu)
+
# Creazione di HDUList e scrittura del file FITS
- hdul = pyfits.HDUList([primary_hdu, ebounds_bin_table_hdu, matrix_bin_table_hdu])
+ hdul = pyfits.HDUList([primary_hdu, matrix_bin_table_hdu, ebounds_bin_table_hdu])
# ! *************************************************************************************************************
diff --git a/python/zero.CANCELLA.py b/python/zero.CANCELLA.py
new file mode 100644
index 0000000000000000000000000000000000000000..abd69086dae30b04a6a88310c486a1aa0efa50ed
--- /dev/null
+++ b/python/zero.CANCELLA.py
@@ -0,0 +1,16 @@
+# Nome del file da leggere
+file_path = "mask.dat"
+
+# Apri il file in modalità lettura
+with open(file_path, 'r') as file:
+ # Itera su ogni riga del file
+ for line_num, line in enumerate(file, start=1):
+ # Rimuovi eventuali spazi vuoti ai bordi e dividi la riga in elementi
+ elements = line.strip().split()
+
+ # Conta il numero di '0' nella riga corrente
+ zero_count = elements.count('0')
+
+ # Stampa il numero di zeri per la riga corrente
+ print(f"Riga {line_num}: {zero_count} zeri")
+
diff --git a/src/DetectorConstruction.cc b/src/DetectorConstruction.cc
index b490a5262764a372216908abe9e8dc678ea0a5c0..7303de24a7e22ff7bea2c10ce5df63ef380839c6 100644
--- a/src/DetectorConstruction.cc
+++ b/src/DetectorConstruction.cc
@@ -39,6 +39,7 @@
#include "G4GDMLParser.hh"
+#include <iostream>
#include <fstream>
#include <vector>
#include <string>
@@ -241,15 +242,15 @@ void DetectorConstruction::DefineMaterials()
G4Element* Cu = new G4Element("Copper" , "Cu", z = 29., a = 63.55*g/mole);
G4Element* Ge = new G4Element("Germanium" , "Ge", z = 32., a = 72.63*g/mole);
// G4Element* As = new G4Element("Arsenic" , "As", z = 33., a = 74.92*g/mole);
- G4Element* Br = new G4Element("Bromine" , "Br", z = 35., a = 79.90*g/mole);
+ G4Element* Br = new G4Element("Bromine" , "Br", z = 35., a = 79.90*g/mole);
// G4Element* Mo = new G4Element("Molybdenum", "Mo", z = 42., a = 95.94*g/mole);
G4Element* I = new G4Element("Iodine" , "I" , z = 53., a = 126.90*g/mole);
G4Element* Cs = new G4Element("Caesium" , "Cs", z = 55., a = 132.90*g/mole);
G4Element* La = new G4Element("Lanthanium", "La", z = 57., a = 138.90*g/mole);
// G4Element* Ta = new G4Element("Tantalum" , "Ta", z = 73., a = 180.94*g/mole);
- G4Element* W = new G4Element("Tungsten" , "W" , z = 74., a = 183.84*g/mole);
+ G4Element* W = new G4Element("Tungsten" , "W" , z = 74., a = 183.84*g/mole);
// G4Element* Pb = new G4Element("Lead" , "Pb", z = 82., a = 207.20*g/mole);
- G4Element* Bi = new G4Element("Bismuth" , "Bi", z = 83., a = 208.98*g/mole);
+ G4Element* Bi = new G4Element("Bismuth" , "Bi", z = 83., a = 208.98*g/mole);
// Materials
// Vacuum
@@ -310,26 +311,26 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
G4SolidStore::GetInstance()->Clean();
- // *** Experimental hall (world volume) ***
+ // ! Experimental hall (world volume)
// Solid
G4Box* experimentalHall_box = new G4Box("experimentalHall_box",
- worldSide/2., worldSide/2., worldSide/2.);
+ worldSide/2., worldSide/2., worldSide/2.);
// Logical
experimentalHall_log = new G4LogicalVolume(experimentalHall_box,
- defaultMaterial, "experimentalHall_log", 0, 0, 0);
+ defaultMaterial, "experimentalHall_log", 0, 0, 0);
// Physical
experimentalHall_phys = new G4PVPlacement(0,
- G4ThreeVector(0,0,0),
- experimentalHall_log,
- "experimentalHall_phys",
- 0,
- false,
- 0);
+ G4ThreeVector(0,0,0),
+ experimentalHall_log,
+ "experimentalHall_phys",
+ 0,
+ false,
+ 0);
- // *** Scintillator ***
+ // ! Scintillator
// Module parameterisation
G4int numberOfBars = n_side*n_side;
@@ -340,9 +341,9 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
G4double scint_container_side_y = (scint_y + scint_gap)*n_side;
G4double scint_container_side_z = scint_z;
G4Box* scint_container_box = new G4Box("scint_container_box",
- scint_container_side_x/2., scint_container_side_y/2., scint_container_side_z/2.);
+ scint_container_side_x/2., scint_container_side_y/2., scint_container_side_z/2.);
G4LogicalVolume* scint_container_log = new G4LogicalVolume(scint_container_box,
- defaultMaterial, "scint_container_log", 0, 0, 0);
+ defaultMaterial, "scint_container_log", 0, 0, 0);
G4VPhysicalVolume* scint_container_phys = new G4PVPlacement(0,
G4ThreeVector(0.,0.,-scint_z/2.),
scint_container_log,
@@ -368,7 +369,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
- // *** Side Collimator ***
+ // ! Side Collimator
G4double coll_side_x = scint_container_side_x;
G4double coll_side_y = scint_container_side_y;
@@ -382,12 +383,12 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
G4double coll_displacement_y = 0;
G4double coll_displacement_z = -scint_z + coll_side_height/2.;
coll_sideXp_phys = new G4PVPlacement(0,
- G4ThreeVector(coll_displacement_x,coll_displacement_y,coll_displacement_z),
- coll_sideXp_log,
- "coll_sideXp_phys",
- experimentalHall_log,
- false,
- 0);
+ G4ThreeVector(coll_displacement_x,coll_displacement_y,coll_displacement_z),
+ coll_sideXp_log,
+ "coll_sideXp_phys",
+ experimentalHall_log,
+ false,
+ 0);
// Solid
G4Box* coll_sideXm_box = new G4Box("coll_sideXm_box", coll_side_thick/2., coll_side_x/2., coll_side_height/2.);
@@ -398,12 +399,12 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
coll_displacement_y = 0.;
coll_displacement_z = -scint_z + coll_side_height/2.;
coll_sideXm_phys = new G4PVPlacement(0,
- G4ThreeVector(coll_displacement_x,coll_displacement_y,coll_displacement_z),
- coll_sideXm_log,
- "coll_sideXm_phys",
- experimentalHall_log,
- false,
- 0);
+ G4ThreeVector(coll_displacement_x,coll_displacement_y,coll_displacement_z),
+ coll_sideXm_log,
+ "coll_sideXm_phys",
+ experimentalHall_log,
+ false,
+ 0);
// Solid
@@ -415,12 +416,12 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
coll_displacement_y = (coll_side_y/2. + coll_side_thick/2.);
coll_displacement_z = -scint_z + coll_side_height/2.;
coll_sideYp_phys = new G4PVPlacement(0,
- G4ThreeVector(coll_displacement_x,coll_displacement_y,coll_displacement_z),
- coll_sideYp_log,
- "coll_sideYp_phys",
- experimentalHall_log,
- false,
- 0);
+ G4ThreeVector(coll_displacement_x,coll_displacement_y,coll_displacement_z),
+ coll_sideYp_log,
+ "coll_sideYp_phys",
+ experimentalHall_log,
+ false,
+ 0);
// Solid
G4Box* coll_sideYm_box = new G4Box("coll_sideYm_box", coll_side_y/2.+coll_side_thick, coll_side_thick/2., coll_side_height/2.);
@@ -431,20 +432,20 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
coll_displacement_y = -(coll_side_y/2. + coll_side_thick/2.);
coll_displacement_z = -scint_z + coll_side_height/2.;
coll_sideYm_phys = new G4PVPlacement(0,
- G4ThreeVector(coll_displacement_x,coll_displacement_y,coll_displacement_z),
- coll_sideYm_log,
- "coll_sideYm_phys",
- experimentalHall_log,
- false,
- 0);
+ G4ThreeVector(coll_displacement_x,coll_displacement_y,coll_displacement_z),
+ coll_sideYm_log,
+ "coll_sideYm_phys",
+ experimentalHall_log,
+ false,
+ 0);
- // *** Spacecraft bus (box approximation) ***
+ // ! Spacecraft bus (box approximation)
// TODO: insert here new S/C model (hexagonal, rotated wrt XGIS unit)
// Solid
G4Box* bus1_box = new G4Box("bus1_box",
- (coll_side_x + 2*coll_side_thick)/2., (coll_side_y + 2*coll_side_thick)/2., bus_thick/2.);
+ (coll_side_x + 2*coll_side_thick)/2., (coll_side_y + 2*coll_side_thick)/2., bus_thick/2.);
// Logical
bus1_log = new G4LogicalVolume(bus1_box, busMaterial, "bus1_log", 0, 0, 0);
// Physical
@@ -452,12 +453,12 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
coll_displacement_y = 0;
coll_displacement_z = -scint_z - optCoupler_thick - sdd_thick - bus_thick/2.;
bus1_phys = new G4PVPlacement(0,
- G4ThreeVector(coll_displacement_x,coll_displacement_y,coll_displacement_z),
- bus1_log,
- "bus1_phys",
- experimentalHall_log,
- false,
- 0);
+ G4ThreeVector(coll_displacement_x,coll_displacement_y,coll_displacement_z),
+ bus1_log,
+ "bus1_phys",
+ experimentalHall_log,
+ false,
+ 0);
// // Solid
// G4Box* bus1_box = new G4Box("bus1_box",
@@ -478,7 +479,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
// 0);
- // *** SDDs ***
+ // ! SDDs
// Solid
G4Box* sdd_box = new G4Box("sdd_box", coll_side_x/2., coll_side_y/2., sdd_thick/2.);
@@ -491,26 +492,26 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
G4double sddUpper_y = 0;
G4double sddUpper_z = sdd_thick/2. + optCoupler_thick;
sddUpper_phys = new G4PVPlacement(0,
- G4ThreeVector(sddUpper_x,sddUpper_y,sddUpper_z),
- sddUpper_log,
- "sddUpper_phys",
- experimentalHall_log,
- false,
- 0);
+ G4ThreeVector(sddUpper_x,sddUpper_y,sddUpper_z),
+ sddUpper_log,
+ "sddUpper_phys",
+ experimentalHall_log,
+ false,
+ 0);
G4double sddLower_x = 0;
G4double sddLower_y = 0;
G4double sddLower_z = -scint_z - sdd_thick/2.;
sddLower_phys = new G4PVPlacement(0,
- G4ThreeVector(sddLower_x,sddLower_y,sddLower_z),
- sddLower_log,
- "sddLower_phys",
- experimentalHall_log,
- false,
- 0);
+ G4ThreeVector(sddLower_x,sddLower_y,sddLower_z),
+ sddLower_log,
+ "sddLower_phys",
+ experimentalHall_log,
+ false,
+ 0);
- // *** Optical coupler ***
+ // ! Optical coupler
// Solid
G4Box* optCoupler_box = new G4Box("sdd_box", coll_side_x/2., coll_side_y/2., sdd_thick/2.);
@@ -522,54 +523,89 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
G4double optCoupler_y = 0;
G4double optCoupler_z = optCoupler_thick/2.;
optCoupler_phys = new G4PVPlacement(0,
- G4ThreeVector(optCoupler_x,optCoupler_y,optCoupler_z),
- optCoupler_log,
- "optCoupler_phys",
- experimentalHall_log,
- false,
- 0);
+ G4ThreeVector(optCoupler_x,optCoupler_y,optCoupler_z),
+ optCoupler_log,
+ "optCoupler_phys",
+ experimentalHall_log,
+ false,
+ 0);
- // *** Mask ***
+ // ! Mask
+
// Solid
G4Box* mask_container_box = new G4Box("mask_box", mask_side_x/2., mask_side_y/2., mask_thick/2.);
+
// Logical
G4LogicalVolume* mask_container_log = new G4LogicalVolume(mask_container_box, defaultMaterial, "mask_container_log", 0, 0, 0);
// Physical
G4VPhysicalVolume* mask_container_phys = new G4PVPlacement(0,
- G4ThreeVector(0,0,mask_focalDistance+sdd_thick+optCoupler_thick),
- mask_container_log,
- "mask_container_phys",
- experimentalHall_log,
- false,
- 0);
+ G4ThreeVector(0,0,mask_focalDistance+sdd_thick+optCoupler_thick),
+ mask_container_log,
+ "mask_container_phys",
+ experimentalHall_log,
+ false,
+ 0);
- G4double mask_elem_side_x = mask_side_x/(2*n_mask_elem_side);
- G4double mask_elem_side_y = mask_side_y/(2*n_mask_elem_side);
-
- G4int n_mask_elem_replicas = (2*n_mask_elem_side) * (2*n_mask_elem_side) * 1/2;
+ // G4double mask_elem_side_x = mask_side_x/(2*n_mask_elem_side);
+ // G4double mask_elem_side_y = mask_side_y/(2*n_mask_elem_side);
+ // G4int n_mask_elem_replicas = (2*n_mask_elem_side) * (2*n_mask_elem_side) * 1/2;
+
+ G4double mask_elem_side_x = 1.027*cm;
+ G4double mask_elem_side_y = 1.027*cm;
+ //G4double mask_elem_side_x = (mask_side_x/n_mask_elem_side)*cm;
+ //G4double mask_elem_side_y = (mask_side_y/n_mask_elem_side)*cm;
+ G4double gap = (mask_side_x/n_mask_elem_side) - mask_elem_side_x;
// Solid
G4Box* mask_box = new G4Box("mask_box", mask_elem_side_x/2., mask_elem_side_y/2., mask_thick/2.);
+
// Logical
mask_log = new G4LogicalVolume(mask_box, maskMaterial, "mask_log", 0, 0, 0);
+
+ // Creazione dell'array per le posizioni degli zeri
+ std::vector<std::pair<int, int>> zeroPositions;
+
// Physical
- maskParam = new MaskParameterisation(mask_elem_side_x, mask_elem_side_y, n_mask_elem_side);
+ MaskParameterisation* maskParam = new MaskParameterisation(mask_elem_side_x, mask_elem_side_y, n_mask_elem_side, gap, "mask.dat");
+
+ int totalZeros = maskParam->GetTotalZeros();
+
+ if (totalZeros != -1) {
+ // for (const auto& pos : zeroPositions) {
+ // G4cout << "Zero a posizione: (" << pos.first << ", " << pos.second << ")" << G4endl;
+ // }
+ // G4cout << "Numero totale di zeri: " << totalZeros << G4endl;
+ // G4cout << "gap: " << gap << G4endl;
+
+ // Posizionamento parametrizzato delle repliche
+ mask_phys = new G4PVParameterised(
+ "mask_phys", // Nome
+ mask_log, // Logica dell'elemento maschera
+ mask_container_log, // Volume madre
+ kXAxis, // Asse lungo cui avviene il posizionamento
+ totalZeros, // Numero di repliche
+ maskParam // Parametrizzazione
+ );
+ } else {
+ G4cerr << "Errore durante il caricamento del file mask.dat." << G4endl;
+ }
+
- mask_phys = new G4PVParameterised("mask_phys", // their name
- mask_log, // their logical volume
- mask_container_log, // Mother logical volume
- kXAxis, // Are placed along this axis
- n_mask_elem_replicas, // Number of replicas
- maskParam); // The parametrisation
+ // mask_phys = new G4PVParameterised("mask_phys", // their name
+ // mask_log, // their logical volume
+ // mask_container_log, // Mother logical volume
+ // kXAxis, // Are placed along this axis
+ // n_mask_elem_replicas, // Number of replicas
+ // maskParam); // The parametrisation
- // *** Upper Collimator, internal ***
+ // ! Upper Collimator, internal
// Slabs with trapezoidal shape, around the SDD
// Remember the calling sequence for the G4Trap class:
@@ -593,14 +629,14 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
G4double base_side_x = coll_side_x + 2*coll_side_thick;
G4double base_side_y = coll_side_y + 2*coll_side_thick;
- G4cout << "base_side_x, base_side_y (cm) " << base_side_x/cm << " " << base_side_y/cm << G4endl;
+ //G4cout << "base_side_x, base_side_y (cm) " << base_side_x/cm << " " << base_side_y/cm << G4endl;
// Rotation angles for the y-side and x-side collimators (i.e. inclination)
double xrotTheta = atan((mask_side_x + coll_upper1_thick/2 - base_side_x)/(2*mask_focalDistance));
double yrotTheta = atan((mask_side_y + coll_upper1_thick/2 - base_side_y)/(2*mask_focalDistance));
- G4cout << "xrotTheta, yrotTheta (rad) " << xrotTheta << " " << yrotTheta << G4endl;
- G4cout << "xrotTheta, yrotTheta (deg) " << 180/3.141592*(xrotTheta) << " " << 180/3.141592*(yrotTheta) << G4endl;
+ //G4cout << "xrotTheta, yrotTheta (rad) " << xrotTheta << " " << yrotTheta << G4endl;
+ //G4cout << "xrotTheta, yrotTheta (deg) " << 180/3.141592*(xrotTheta) << " " << 180/3.141592*(yrotTheta) << G4endl;
// Solid (for yPlus, yMinus)
G4double pDzY = mask_focalDistance/cos(yrotTheta)/2. ; // half-height
@@ -611,7 +647,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
G4double pDx3 = (base_side_x + (mask_side_x - base_side_x))/2.; // mask side
G4double pDx4 = pDx3; // the other side of the upper face is the same
G4Trap* coll_upper1_trap_y = new G4Trap("coll_upper1_trap_y",
- pDzY, 0, 0, pDy1, pDx1, pDx2, 0, pDy2, pDx3, pDx4, 0);
+ pDzY, 0, 0, pDy1, pDx1, pDx2, 0, pDy2, pDx3, pDx4, 0);
// Solid (for xPlus, xMinus)
G4double pDzX = mask_focalDistance/cos(xrotTheta)/2. ; // half-height
@@ -620,17 +656,17 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
pDx3 = (base_side_y + (mask_side_y - base_side_y))/2.; // mask side
pDx4 = pDx3;
G4Trap* coll_upper1_trap_x = new G4Trap("coll_upper1_trap_x",
- pDzX, 0, 0, pDy1, pDx1, pDx2, 0, pDy2, pDx3, pDx4, 0);
+ pDzX, 0, 0, pDy1, pDx1, pDx2, 0, pDy2, pDx3, pDx4, 0);
// Logical
coll_upper1Yp_log = new G4LogicalVolume(coll_upper1_trap_y,
- collUpper1Material, "coll_upper1Yp_log", 0, 0, 0);
+ collUpper1Material, "coll_upper1Yp_log", 0, 0, 0);
coll_upper1Ym_log = new G4LogicalVolume(coll_upper1_trap_y,
- collUpper1Material, "coll_upper1Ym_log", 0, 0, 0);
+ collUpper1Material, "coll_upper1Ym_log", 0, 0, 0);
coll_upper1Xp_log = new G4LogicalVolume(coll_upper1_trap_x,
- collUpper1Material, "coll_upper1Xp_log", 0, 0, 0);
+ collUpper1Material, "coll_upper1Xp_log", 0, 0, 0);
coll_upper1Xm_log = new G4LogicalVolume(coll_upper1_trap_x,
- collUpper1Material, "coll_upper1Xm_log", 0, 0, 0);
+ collUpper1Material, "coll_upper1Xm_log", 0, 0, 0);
// Physical
// yPlus collimator
@@ -638,24 +674,24 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
yRotPlus->rotateX(-yrotTheta);
G4ThreeVector yShiftPlus(0, base_side_y/2.+ coll_upper1_thick/2. + pDzY*sin(yrotTheta), pDzY*cos(yrotTheta));
coll_upper1Yp_phys = new G4PVPlacement(
- G4Transform3D(*yRotPlus, yShiftPlus),
- coll_upper1Yp_log,
- "collUpper1Yp_phys",
- experimentalHall_log,
- false,
- 0);
+ G4Transform3D(*yRotPlus, yShiftPlus),
+ coll_upper1Yp_log,
+ "collUpper1Yp_phys",
+ experimentalHall_log,
+ false,
+ 0);
// yMinus collimator
G4RotationMatrix* yRotMinus = new G4RotationMatrix;
yRotMinus->rotateX(yrotTheta);
G4ThreeVector yShiftMinus(0, -(base_side_y/2.+ coll_upper1_thick/2. + pDzY*sin(yrotTheta)), pDzY*cos(yrotTheta));
coll_upper1Ym_phys = new G4PVPlacement(
- G4Transform3D(*yRotMinus, yShiftMinus),
- coll_upper1Ym_log,
- "collUpper1Ym_phys",
- experimentalHall_log,
- false,
- 0);
+ G4Transform3D(*yRotMinus, yShiftMinus),
+ coll_upper1Ym_log,
+ "collUpper1Ym_phys",
+ experimentalHall_log,
+ false,
+ 0);
// xPlus collimator
G4RotationMatrix* xRotPlus = new G4RotationMatrix;
@@ -663,12 +699,12 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
xRotPlus->rotateY(xrotTheta);
G4ThreeVector xShiftPlus(base_side_x/2.+ coll_upper1_thick/2. + pDzX*sin(xrotTheta), 0, pDzX*cos(xrotTheta));
coll_upper1Xp_phys = new G4PVPlacement(
- G4Transform3D(*xRotPlus, xShiftPlus),
- coll_upper1Xp_log,
- "collUpper1Xp_phys",
- experimentalHall_log,
- false,
- 0);
+ G4Transform3D(*xRotPlus, xShiftPlus),
+ coll_upper1Xp_log,
+ "collUpper1Xp_phys",
+ experimentalHall_log,
+ false,
+ 0);
// xMinus collimator
G4RotationMatrix* xRotMinus = new G4RotationMatrix;
@@ -676,15 +712,15 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
xRotMinus->rotateY(-xrotTheta);
G4ThreeVector xShiftMinus(-(base_side_x/2.+ coll_upper1_thick/2. + pDzX*sin(xrotTheta)), 0, pDzX*cos(xrotTheta));
coll_upper1Xm_phys = new G4PVPlacement(
- G4Transform3D(*xRotMinus, xShiftMinus),
- coll_upper1Xm_log,
- "collUpper1Xm_phys",
- experimentalHall_log,
- false,
- 0);
+ G4Transform3D(*xRotMinus, xShiftMinus),
+ coll_upper1Xm_log,
+ "collUpper1Xm_phys",
+ experimentalHall_log,
+ false,
+ 0);
- // *** Upper Collimator, external ***
+ // ! Upper Collimator, external
// Slabs with trapezoidal shape, around the SDD
// Remember the calling sequence for the G4Trap class:
@@ -721,7 +757,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
G4double p2Dx3 = (base_side_x + (scaled_mask_side_x - base_side_x))/2.; // mask side
G4double p2Dx4 = p2Dx3; // the other side of the upper face is the same
G4Trap* coll_upper2_trap_y = new G4Trap("coll_upper2_trap_y",
- p2DzY, 0, 0, p2Dy1, p2Dx1, p2Dx2, 0, p2Dy2, p2Dx3, p2Dx4, 0);
+ p2DzY, 0, 0, p2Dy1, p2Dx1, p2Dx2, 0, p2Dy2, p2Dx3, p2Dx4, 0);
// Solid (for xPlus, xMinus)
G4double p2DzX = coll_upper2_height/cos(xrotTheta)/2. ; // half-height
@@ -730,7 +766,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
p2Dx3 = (base_side_y + (scaled_mask_side_y - base_side_y))/2.; // mask side
p2Dx4 = p2Dx3;
G4Trap* coll_upper2_trap_x = new G4Trap("coll_upper2_trap_x",
- p2DzX, 0, 0, p2Dy1, p2Dx1, p2Dx2, 0, p2Dy2, p2Dx3, p2Dx4, 0);
+ p2DzX, 0, 0, p2Dy1, p2Dx1, p2Dx2, 0, p2Dy2, p2Dx3, p2Dx4, 0);
// Logical
coll_upper2Yp_log = new G4LogicalVolume(coll_upper2_trap_y,
@@ -746,42 +782,42 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
// yPlus collimator
G4ThreeVector yShiftPlus2(0, base_side_y/2.+ coll_upper1_thick/cos(yrotTheta) + coll_upper2_thick/2. + p2DzY*sin(yrotTheta), p2DzY*cos(yrotTheta));
coll_upper2Yp_phys = new G4PVPlacement(
- G4Transform3D(*yRotPlus, yShiftPlus2),
- coll_upper2Yp_log,
- "collUpper2Yp_phys",
- experimentalHall_log,
- false,
- 0);
+ G4Transform3D(*yRotPlus, yShiftPlus2),
+ coll_upper2Yp_log,
+ "collUpper2Yp_phys",
+ experimentalHall_log,
+ false,
+ 0);
// yMinus collimator
G4ThreeVector yShiftMinus2(0, -(base_side_y/2.+ coll_upper1_thick/cos(yrotTheta) + coll_upper2_thick/2. + p2DzY*sin(yrotTheta)), p2DzY*cos(yrotTheta));
coll_upper2Ym_phys = new G4PVPlacement(
- G4Transform3D(*yRotMinus, yShiftMinus2),
- coll_upper2Ym_log,
- "collUpper2Ym_phys",
- experimentalHall_log,
- false,
- 0);
+ G4Transform3D(*yRotMinus, yShiftMinus2),
+ coll_upper2Ym_log,
+ "collUpper2Ym_phys",
+ experimentalHall_log,
+ false,
+ 0);
// xPlus collimator
G4ThreeVector xShiftPlus2(base_side_x/2.+ coll_upper1_thick/cos(xrotTheta) + coll_upper2_thick/2. + p2DzX*sin(xrotTheta), 0, p2DzX*cos(xrotTheta));
coll_upper2Xp_phys = new G4PVPlacement(
- G4Transform3D(*xRotPlus, xShiftPlus2),
- coll_upper2Xp_log,
- "collUpper2Xp_phys",
- experimentalHall_log,
- false,
- 0);
+ G4Transform3D(*xRotPlus, xShiftPlus2),
+ coll_upper2Xp_log,
+ "collUpper2Xp_phys",
+ experimentalHall_log,
+ false,
+ 0);
// xMinus collimator
G4ThreeVector xShiftMinus2(-(base_side_x/2.+ coll_upper1_thick/cos(xrotTheta) + coll_upper2_thick/2. + p2DzX*sin(xrotTheta)), 0, p2DzX*cos(xrotTheta));
coll_upper2Xm_phys = new G4PVPlacement(
- G4Transform3D(*xRotMinus, xShiftMinus2),
- coll_upper2Xm_log,
- "collUpper2Xm_phys",
- experimentalHall_log,
- false,
- 0);
+ G4Transform3D(*xRotMinus, xShiftMinus2),
+ coll_upper2Xm_log,
+ "collUpper2Xm_phys",
+ experimentalHall_log,
+ false,
+ 0);
@@ -844,14 +880,15 @@ G4VPhysicalVolume* DetectorConstruction::Construct()
// Dump the geometry in a GDML file, solo se non esiste
G4String fWriteFile = "geometry.gdml";
- std::ifstream fileCheck(fWriteFile.c_str());
- if (fileCheck.good()) {
- std::cout << "il file geometry.gdml esiste già" << std::endl;
- } else {
- G4GDMLParser parser;
- //std::remove(fWriteFile.c_str());
- parser.Write(fWriteFile, experimentalHall_phys, false);
- }
+ // questo pezzo non va bene perche se modifico non me lo sovrascrivr //! MIGLIORARE
+ // std::ifstream fileCheck(fWriteFile.c_str());
+ // if (fileCheck.good()) {
+ // G4cout << "il file geometry.gdml esiste già" << G4endl;
+ // } else {
+ // G4GDMLParser parser;
+ // //std::remove(fWriteFile.c_str());
+ // parser.Write(fWriteFile, experimentalHall_phys, false);
+ // }
// The function must return the physical volume of the world
return experimentalHall_phys;
@@ -873,9 +910,9 @@ void DetectorConstruction::ConstructSDandField()
// Instantiation of the sensitive detector and readout geometry
SensitiveDetector* scint_SD = new SensitiveDetector("SCI",
- n_side,
- scint_x+scint_gap, scint_z,
- NbOfSlicesPerBar);
+ n_side,
+ scint_x+scint_gap, scint_z,
+ NbOfSlicesPerBar);
sdman->AddNewDetector(scint_SD); // Mandatory since Geant v. 4.10.03
@@ -883,8 +920,8 @@ void DetectorConstruction::ConstructSDandField()
// Instantiation of the sensitive detector and readout geometry
SDDSensitiveDetector* sdd_SD = new SDDSensitiveDetector("SDD",
- scint_x+scint_gap,
- n_side, n_side);
+ scint_x+scint_gap,
+ n_side, n_side);
sdman->AddNewDetector(sdd_SD); // Mandatory since Geant v. 4.10.03
SetSensitiveDetector(sddUpper_log, sdd_SD);
}
@@ -947,6 +984,9 @@ void DetectorConstruction::PrintParameters()
<< "\n Scint Side X (cm) : " << scint_x/cm
<< "\n Scint Side Y (cm) : " << scint_y/cm
<< "\n Scint Side Z (cm) : " << scint_z/cm
+ << "\n mask_side_x (cm) : " << mask_side_x/cm
+ << "\n mask_side_y (cm) : " << mask_side_y/cm
+ << "\n n_mask_elem_side : " << n_mask_elem_side
<< "\n Scint material : " << scintMaterial -> GetName()
<< "\n Mask material : " << maskMaterial -> GetName()
<< "\n Collimator 1 material : " << collUpper1Material -> GetName()
diff --git a/src/MaskParameterisation.cc b/src/MaskParameterisation.cc
index 0009ab80ba0850f3e811e02779329392c4977603..38f85a220f0a2486a15ce019e2ac34d67bf00b78 100644
--- a/src/MaskParameterisation.cc
+++ b/src/MaskParameterisation.cc
@@ -7,12 +7,21 @@
#include "G4RotationMatrix.hh"
#include "G4Box.hh"
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <sstream>
+#include <vector>
-MaskParameterisation::MaskParameterisation(G4double side_length_X, G4double side_length_Y, G4int n_side)
+
+MaskParameterisation::MaskParameterisation(G4double side_length_X, G4double side_length_Y, G4int n_side, G4double gap, const std::string& filename)
+ : fSideX(side_length_X), fSideY(side_length_Y), nSide(n_side), fGap(gap)
{
- fSideX = side_length_X;
- fSideY = side_length_Y;
- nSide = n_side;
+ // Carica zeroPositions dal file
+ int totalZeros = LoadZeroPositions(filename);
+ if (totalZeros == -1) {
+ G4cerr << "Errore durante il caricamento del file " << filename << G4endl;
+ }
}
@@ -20,26 +29,89 @@ MaskParameterisation::~MaskParameterisation()
{;}
+int MaskParameterisation::LoadZeroPositions(const std::string& filename)
+{
+ std::ifstream infile(filename); //apre il file in modalità lettura
+ if (!infile) {
+ std::cerr << "Errore nell'apertura del file " << filename << std::endl;
+ return -1;
+ }
+
+ std::string line; //memorizza la riga corrente del file
+ int rowID = 0;
+ int totalZeros = 0;
+
+ // Leggi il file riga per riga
+ while (std::getline(infile, line)) { // legge il file una riga alla volta. std::getline carica il contenuto della riga corrente in line
+ // Salta le righe commentate
+ if (line[0] == '#') continue;
+
+ std::istringstream iss(line); // crea uno stream di input (iss) dalla riga letta (line), per processare ogni elemento della riga come valore numerico
+ int colID = 0;
+ int value; // value memorizza ciascun numero estratto dalla riga.
+
+ // Leggi ogni numero della riga e controlla se è uno zero
+ while (iss >> value) { // >> estrae il prossimo elemento dallo stream iss e lo assegna a value.
+ if (value == 0) { // Se l'operazione ha successo (ossia se iss ha ancora un numero da leggere), il ciclo while continua; altrimenti, iss >> value restituisce false, terminando il ciclo while
+ // Aggiungi la posizione (rowID, colID) all'array zeroPositions
+ zeroPositions.push_back(std::make_pair(rowID, colID));
+ totalZeros++;
+ }
+ colID++;
+ }
+ rowID++;
+ }
+
+ infile.close();
+ return totalZeros;
+}
+
+int MaskParameterisation::GetTotalZeros() const
+{
+ return zeroPositions.size(); // Ritorna il numero totale di zeri, che corrisponde alla dimensione del vettore zeroPositions
+}
+
+
void MaskParameterisation::ComputeTransformation (const G4int copyNo, G4VPhysicalVolume* physVol) const
{
G4double xCell = 0.*cm;
G4double yCell = 0.*cm;
G4double zCell = 0.*cm;
- if (nSide % 2 == 0) // Even number of tiles per side
- {
- G4int rowID = copyNo / nSide; // identifies the row of the tile
- G4int colID = (copyNo % nSide)*2 + (((int) copyNo/nSide)%2); // in order to have tiles in odd-numbered rows shifted by one column w.r.t. the adjacent rows
-// G4cout << "*** DEBUG copyNo " << copyNo << " rowID " << rowID << " colID " << colID << G4endl;
- xCell = (2*colID - (2*nSide -1)) *fSideX/2;
- yCell = (2*rowID - (2*nSide -1)) *fSideY/2;
-// G4cout << "*** DEBUG xCell " << xCell/cm << " yCell " << yCell/cm << G4endl;
-// G4cout << " " << G4endl;
- }
-
- physVol->SetTranslation(G4ThreeVector(xCell, yCell, zCell));
+ int rowID = zeroPositions[copyNo].first;
+ int colID = zeroPositions[copyNo].second;
+
+ xCell = (fSideX) *((-nSide+1)/2 + colID); // potrei scrivere fSideX+fGap se voglio aggiungere il gap
+ yCell = (fSideY) *((-nSide+1)/2 + rowID);
+
+ physVol->SetTranslation(G4ThreeVector(xCell, yCell, zCell));
+
+ // G4cout << "*** DEBUG copyNo " << copyNo << " rowID " << rowID << " colID " << colID << G4endl;
+ // G4cout << "*** DEBUG xCell " << xCell/cm << " yCell " << yCell/cm << G4endl;
}
+
+// ! OLD MASK
+// void MaskParameterisation::ComputeTransformation (const G4int copyNo, G4VPhysicalVolume* physVol) const
+// {
+// G4double xCell = 0.*cm;
+// G4double yCell = 0.*cm;
+// G4double zCell = 0.*cm;
+
+// if (nSide % 2 == 0) // Even number of tiles per side
+// {
+// G4int rowID = copyNo / nSide; // identifies the row of the tile
+// G4int colID = (copyNo % nSide)*2 + (((int) copyNo/nSide)%2); // in order to have tiles in odd-numbered rows shifted by one column w.r.t. the adjacent rows
+// // G4cout << "*** DEBUG copyNo " << copyNo << " rowID " << rowID << " colID " << colID << G4endl;
+// xCell = (2*colID - (2*nSide -1)) *fSideX/2;
+// yCell = (2*rowID - (2*nSide -1)) *fSideY/2;
+// // G4cout << "*** DEBUG xCell " << xCell/cm << " yCell " << yCell/cm << G4endl;
+// // G4cout << " " << G4endl;
+// }
+
+// physVol->SetTranslation(G4ThreeVector(xCell, yCell, zCell));
+
+// }
diff --git a/xgis_M7 b/xgis_M7
index c03587f07ce55c209a5acd883d4afa7742510003..7c7706bbcc79ec7141751625e5ddc9cc7d80623c 100755
Binary files a/xgis_M7 and b/xgis_M7 differ