Skip to content

Argon #257

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
kvjmistry wants to merge 14 commits into
base: master
Choose a base branch
from
Open

Argon #257

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
18ede98
Update argon density calculation
kvjmistry Jun 7, 2024
d29cbe2
Add argon gas to NEXT100 geometry
kvjmistry Jun 7, 2024
6aa3b7e
Add configurable drift vel
kvjmistry Jun 13, 2024
37fe9f3
Update Ar gas density interpolation
kvjmistry Jun 13, 2024
4faaae5
Add example full sim alpha config
kvjmistry Jun 13, 2024
a3e2475
Cosmetics
kvjmistry Jun 13, 2024
dce2a8b
Fix typo
kvjmistry Jun 13, 2024
2444b95
Revert argon definition
kvjmistry Jun 17, 2024
61b56e9
Add new mm per mus unit to drift vel
kvjmistry Jun 17, 2024
c468342
Cosmetics
kvjmistry Jun 19, 2024
5bf2f31
Update refractive index for FS and Sapphire for high E VUV photons
kvjmistry Apr 7, 2025
f1b86b2
Update N100 Optical Geom for use with Argon
kvjmistry Apr 7, 2025
6b9d890
Fix typos in verbosity
kvjmistry Apr 7, 2025
9984802
Update alpha config
kvjmistry Apr 7, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
64 changes: 64 additions & 0 deletions macros/NEXT100_Ar_alpha_full.config.mac
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
## ----------------------------------------------------------------------------
## nexus | NEXT100_Ar_alpha_full.config.mac
##
## Initialization macro to alpha particle decays in the active
## volume of the NEXT-100 detector during the low pressure run.
##
## The NEXT Collaboration
## ----------------------------------------------------------------------------

##### VERBOSITY #####
/run/verbose 0
/event/verbose 0
/tracking/verbose 0

/process/em/verbose 0

##### GEOMETRY #####
/Geometry/Next100/elfield true
/Geometry/Next100/EL_field 6.7 kV/cm
/Geometry/Next100/max_step_size 1. mm
/Geometry/Next100/pressure 4.15 bar
/Geometry/Next100/gas GAr
/Geometry/PmtR11410/time_binning 25. nanosecond
/Geometry/Next100/sipm_time_binning 1. microsecond

/Geometry/Next100/drift_v 1.480 mm/microsecond
/Geometry/Next100/EL_drift_v 6.249 mm/microsecond

/Geometry/Next100/drift_transv_diff 1.87 mm/sqrt(cm)
/Geometry/Next100/drift_long_diff 0.65 mm/sqrt(cm)

/Geometry/Next100/ELtransv_diff 0.46 mm/sqrt(cm)
/Geometry/Next100/ELlong_diff 0.30 mm/sqrt(cm)

/Geometry/Next100/e_lifetime 50 ms

##### GENERATOR #####
/Generator/SingleParticle/particle alpha

# For N100 use 5.6 MeV here
/Generator/SingleParticle/min_energy 0.01 MeV
/Generator/SingleParticle/max_energy 0.01 MeV
/Generator/SingleParticle/region ACTIVE

##### ACTIONS #####
/Actions/DefaultEventAction/min_energy 50 keV
/Actions/DefaultEventAction/max_energy 6 MeV

##### PHYSICS #####
## Full simulation
/PhysicsList/Nexus/clustering true
/PhysicsList/Nexus/drift true
/PhysicsList/Nexus/electroluminescence true

##### PERSISTENCY #####
/nexus/persistency/output_file NEXT100_Ar_alpha

/nexus/random_seed 100
/nexus/persistency/start_id 100

## eventType options: bb0nu, bb2nu, background
/nexus/persistency/event_type background

/process/optical/processActivation Cerenkov false
Copy link
Collaborator

@gonzaponte gonzaponte Apr 8, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Missing a newline at the end of file.

27 changes: 27 additions & 0 deletions macros/NEXT100_Ar_alpha_full.init.mac
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
## ----------------------------------------------------------------------------
## nexus | NEXT100_Ar_alpha_full.init.mac
##
## Initialization macro to alpha particle decays in the active
## volume of the NEXT-100 detector during the low pressure run.
##
## The NEXT Collaboration
## ----------------------------------------------------------------------------

/PhysicsList/RegisterPhysics G4EmStandardPhysics_option4
/PhysicsList/RegisterPhysics G4DecayPhysics
/PhysicsList/RegisterPhysics G4RadioactiveDecayPhysics
/PhysicsList/RegisterPhysics NexusPhysics
/PhysicsList/RegisterPhysics G4StepLimiterPhysics
/PhysicsList/RegisterPhysics G4OpticalPhysics

/nexus/RegisterGeometry Next100

/nexus/RegisterGenerator SingleParticleGenerator

/nexus/RegisterPersistencyManager PersistencyManager

/nexus/RegisterRunAction DefaultRunAction
/nexus/RegisterEventAction DefaultEventAction
/nexus/RegisterTrackingAction DefaultTrackingAction

/nexus/RegisterMacro macros/NEXT100_Ar_alpha_full.config.mac
11 changes: 10 additions & 1 deletion source/geometries/Next100FieldCage.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include "OpticalMaterialProperties.h"
#include "UniformElectricDriftField.h"
#include "XenonProperties.h"
#include "ArgonGasProperties.h"
#include "CylinderPointSampler.h"
#include "BoxPointSampler.h"
#include "HexagonMeshTools.h"
Expand Down Expand Up @@ -742,7 +743,15 @@ void Next100FieldCage::BuildELRegion()
el_field->SetDriftVelocity(EL_drift_v_);
el_field->SetTransverseDiffusion(ELtransv_diff_);
el_field->SetLongitudinalDiffusion(ELlong_diff_);
el_field->SetLightYield(XenonELLightYield(ELelectric_field_, pressure_));

// Decide whether to use argon or xenon LY
if (gas_->GetName() == "GAr"){
el_field->SetLightYield(ArgonELLightYield(ELelectric_field_, pressure_));
}
else {
el_field->SetLightYield(XenonELLightYield(ELelectric_field_, pressure_));
}

Comment on lines -745 to +754
Copy link
Collaborator

@gonzaponte gonzaponte Apr 8, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we protect against invalid names like everywhere else?

G4Region* el_region = new G4Region("EL_REGION");
el_region->SetUserInformation(el_field);
el_region->AddRootLogicalVolume(el_gap_logic);
Expand Down
23 changes: 12 additions & 11 deletions source/geometries/Next100OpticalGeometry.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -105,25 +105,26 @@ namespace nexus {

if (gas_ == "naturalXe") {
gas_mat = materials::GXe(pressure_, temperature_);
gas_mat->SetMaterialPropertiesTable(opticalprops::GXe(pressure_,
temperature_,
sc_yield_,
e_lifetime_));
} else if (gas_ == "enrichedXe") {
gas_mat = materials::GXeEnriched(pressure_, temperature_);
gas_mat->SetMaterialPropertiesTable(opticalprops::GXe(pressure_,
temperature_,
sc_yield_,
e_lifetime_));
} else if (gas_ == "depletedXe") {
gas_mat = materials::GXeDepleted(pressure_, temperature_);
} else if (gas_ == "GAr") {
gas_mat = materials::GAr(pressure_, temperature_);
} else {
G4Exception("[Next100OpticalGeometry]", "Construct()", FatalException,
"Unknown kind of gas, valid options are: naturalXe, enrichedXe, depletedXe, GAr.");
}

// Gas Properties
if (gas_ == "GAr"){
gas_mat->SetMaterialPropertiesTable(opticalprops::GAr(sc_yield_, e_lifetime_));
}
else {
gas_mat->SetMaterialPropertiesTable(opticalprops::GXe(pressure_,
temperature_,
sc_yield_,
e_lifetime_));
} else {
G4Exception("[Next100OpticalGeometry]", "Construct()", FatalException,
"Unknown kind of gas, valid options are: naturalXe, enrichedXe, depletedXe.");
}

G4double gas_size = lab_size - 10.*cm;
Expand Down
10 changes: 9 additions & 1 deletion source/geometries/Next100Vessel.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -342,16 +342,24 @@ namespace nexus {
} else if (gas_ == "XeHe") {
vessel_gas_mat = materials::GXeHe(pressure_, 300. * kelvin,
xe_perc_, helium_mass_num_);
} else if (gas_ == "GAr") {
vessel_gas_mat = materials::GAr(pressure_, temperature_);
} else {
G4Exception("[Next100Vessel]", "Construct()", FatalException,
"Unknown kind of xenon, valid options are: "
"natural, enriched, depleted, or XeHe.");
"naturalXe, enrichedXe, depletedXe, GAr, or XeHe.");
}

// Gas Properties
if (gas_ == "GAr"){
vessel_gas_mat->SetMaterialPropertiesTable(opticalprops::GAr(sc_yield_, e_lifetime_));
}
else {
vessel_gas_mat->SetMaterialPropertiesTable(opticalprops::GXe(pressure_,
temperature_,
sc_yield_,
e_lifetime_));
}
Comment on lines +353 to +362
Copy link
Collaborator

@gonzaponte gonzaponte Apr 8, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is an issue for a different PR, but it seems that we are defining the same gas in different places, perhaps we should revisit this and propagate the gas material from the top down.


G4LogicalVolume* vessel_gas_logic =
new G4LogicalVolume(vessel_gas_final_solid, vessel_gas_mat, "VESSEL_GAS");
Expand Down
66 changes: 43 additions & 23 deletions source/materials/ArgonGasProperties.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,29 +24,49 @@ namespace nexus {

G4double ArgonDensity(G4double pressure)
{
//These values are for a temperature of 300 K
// taken from http://www.nist.gov/srd/upload/jpcrd363.pdf
G4double density = 1.60279*kg/m3;

if (pressure/bar > 0.9 && pressure/bar < 1.1)
density = 1.60279*kg/m3;
else if (pressure/bar > 1.9 && pressure/bar < 2.1)
density = 3.20719*kg/m3;
else if (pressure/bar > 4.9 && pressure/bar < 5.1)
density = 8.032*kg/m3;
else if (pressure/bar > 9.9 && pressure/bar < 10.1)
density = 16.1118*kg/m3;
else if (pressure/bar > 14.9 && pressure/bar < 15.1)
density = 24.2369 *kg/m3;
else if (pressure/bar > 19.9 && pressure/bar < 20.1)
density = 32.4066*kg/m3;
else if (pressure/bar > 29.9 && pressure/bar < 30.1)
density = 48.8708*kg/m3;
else if (pressure/bar > 39.9 && pressure/bar < 40.1)
density = 65.494*kg/m3;
else
G4Exception("[ArgonProperties]", "ArgonDensity()", FatalException,
"Unknown argon density for this pressure!");
// Computes Ar (gas) density at T = 293 K
// values taken from https://webbook.nist.gov/chemistry/fluid).
// We assume a linear interpolation between any pair of values in the database.

G4double density;

const G4int n_pressures = 14;
G4double data[n_pressures][2] = {{ 1.0 * bar, 1.641 * kg/m3},
{ 2.0 * bar, 3.284 * kg/m3},
{ 3.0 * bar, 4.929 * kg/m3},
{ 4.0 * bar, 7.402 * kg/m3},
{ 5.0 * bar, 8.2268 * kg/m3},
{ 6.0 * bar, 9.8788 * kg/m3},
{ 7.0 * bar, 11.533 * kg/m3},
{ 8.0 * bar, 13.189 * kg/m3},
{ 9.0 * bar, 14.848 * kg/m3},
{ 10.0 * bar, 16.508 * kg/m3},
{ 15.0 * bar, 24.843 * kg/m3},
{ 20.0 * bar, 33.231 * kg/m3},
{ 30.0 * bar, 50.155 * kg/m3}};
G4bool found = false;
Copy link
Collaborator

@gonzaponte gonzaponte Apr 8, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you could initialize density to 0 and then check if it's still 0 later so you don't need to keep track of a separate variable


for (G4int i=0; i<n_pressures-1; ++i) {
if (pressure >= data[i][0] && pressure < data[i+1][0]) {
G4double x1 = data[i][0];
G4double x2 = data[i+1][0];
G4double y1 = data[i][1];
G4double y2 = data[i+1][1];
density = y1 + (y2-y1)*(pressure-x1)/(x2-x1);
found = true;
break;
}
}

if (!found) {
if (pressure == data[n_pressures-1][0]) {
density = data[n_pressures-1][1];
}
else {
throw "Unknown argon density for this pressure!";
}
}


return density;
}
Expand Down
23 changes: 18 additions & 5 deletions source/materials/OpticalMaterialProperties.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,13 +65,20 @@ namespace opticalprops {
for (int i=0; i<ri_entries; i++) {
rIndex.push_back(seq.RefractiveIndex(hc_/ri_energy[i]));
}
ri_energy.push_back(optPhotMaxE_); // This sets the refractive index between optPhotFusedSilicaMaxE_ and
rIndex.push_back(rIndex[rIndex.size()-1]); // optPhotMaxE_ to the value obtained at optPhotFusedSilicaMaxE_

// The rindex is not defined for sapphire beyond optPhotFusedSilicaMaxE_ = 10.7 eV
// Measurements https://link.springer.com/article/10.1134/S0020441206030195 imply
// that the transmission goes to zero for energies higher than this (>120 nm)
// We set n to be an arbitarily high value of 10 so the value is
// implemented in the simulation (transmission goes to zero)
ri_energy.push_back(optPhotMaxE_);
rIndex.push_back(10);
Comment on lines +69 to +75
Copy link
Collaborator

@gonzaponte gonzaponte Apr 8, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I recall correctly, G4 interpolates between values, so the refractive index between optPhotFusedSilicaMaxE_ and optPhotMaxE_ will grow linearly. Is this what we want or do we want it to be high immediately after optPhotFusedSilicaMaxE_ ?


// for (unsigned int i=0; i<ri_energy.size(); i++) {
// G4cout << "* FusedSilica rIndex: " << std::setw(5) << ri_energy[i]/eV
// << " eV -> " << rIndex[i] << G4endl;
// }

mpt->AddProperty("RINDEX", ri_energy, rIndex);

// ABSORPTION LENGTH
Expand Down Expand Up @@ -431,14 +438,20 @@ namespace opticalprops {
for (int i=0; i<ri_entries; i++) {
rIndex.push_back(seq.RefractiveIndex(hc_/ri_energy[i]));
}
// This sets the refractive index between optPhotSapphireMaxE_ and
// optPhotMaxE_ to the value obtained at optPhotSapphireMaxE_

// The rindex is not defined for sapphire beyond optPhotSapphireMaxE_ = 10.3 eV
// Measurements https://link.springer.com/article/10.1134/S0020441206030195 imply
// that the transmission goes to zero for energies higher than this (>120 nm)
// We set n to be an arbitarily high value of 10 so the value is
// implemented in the simulation (transmission goes to zero)
ri_energy.push_back(optPhotMaxE_);
rIndex.push_back(rIndex[rIndex.size()-1]);
rIndex.push_back(10);

Comment on lines -434 to +449
Copy link
Collaborator

@gonzaponte gonzaponte Apr 8, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same comment as above.

// for (unsigned int i=0; i<ri_energy.size(); i++) {
// G4cout << "* Sapphire rIndex: " << std::setw(5)
// << ri_energy[i]/eV << " eV -> " << rIndex[i] << G4endl;
// }

mpt->AddProperty("RINDEX", ri_energy, rIndex);

// ABSORPTION LENGTH
Expand Down