diff --git a/macros/NEXT100_Ar_alpha_full.config.mac b/macros/NEXT100_Ar_alpha_full.config.mac new file mode 100644 index 000000000..fc61c96db --- /dev/null +++ b/macros/NEXT100_Ar_alpha_full.config.mac @@ -0,0 +1,56 @@ +## ---------------------------------------------------------------------------- +## 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.2 bar +/Geometry/Next100/gas GAr +/Geometry/PmtR11410/time_binning 25. nanosecond +/Geometry/Next100/sipm_time_binning 1. microsecond + +## Geant4 only permits units of mm/ns via the macro +/Geometry/Next100/drift_vel 1.480 mm/mus +/Geometry/Next100/ELdrift_vel 6.249 mm/mus + +/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) + + +##### GENERATOR ##### +/Generator/SingleParticle/particle alpha +/Generator/SingleParticle/min_energy 5.6 MeV +/Generator/SingleParticle/max_energy 5.6 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 +## eventType options: bb0nu, bb2nu, background +/nexus/persistency/event_type background diff --git a/macros/NEXT100_Ar_alpha_full.init.mac b/macros/NEXT100_Ar_alpha_full.init.mac new file mode 100644 index 000000000..b2e5f9d2a --- /dev/null +++ b/macros/NEXT100_Ar_alpha_full.init.mac @@ -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 diff --git a/source/geometries/Next100FieldCage.cc b/source/geometries/Next100FieldCage.cc index 94a334983..7a203328d 100644 --- a/source/geometries/Next100FieldCage.cc +++ b/source/geometries/Next100FieldCage.cc @@ -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" @@ -99,6 +100,9 @@ Next100FieldCage::Next100FieldCage(G4double grid_thickn): drift_long_diff_ (.3 * mm/sqrt(cm)), ELtransv_diff_ (0. * mm/sqrt(cm)), ELlong_diff_ (0. * mm/sqrt(cm)), + // Drift velocities + drift_vel_ (1.0 * mm/microsecond), + ELdrift_vel_ (2.5 * mm/microsecond), // EL electric field elfield_ (0), ELelectric_field_ (34.5*kilovolt/cm), @@ -117,6 +121,7 @@ Next100FieldCage::Next100FieldCage(G4double grid_thickn): /// Define new categories new G4UnitDefinition("kilovolt/cm","kV/cm","Electric field", kilovolt/cm); new G4UnitDefinition("mm/sqrt(cm)","mm/sqrt(cm)","Diffusion", mm/sqrt(cm)); + new G4UnitDefinition("mm/mus","mm/mus","Velocity", mm/microsecond); /// Initializing the geometry navigator (used in vertex generation) geom_navigator_ = @@ -155,6 +160,19 @@ Next100FieldCage::Next100FieldCage(G4double grid_thickn): ELlong_diff_cmd.SetParameterName("ELlong_diff", true); ELlong_diff_cmd.SetUnitCategory("Diffusion"); + G4GenericMessenger::Command& drift_velocity_cmd = + msg_->DeclareProperty("drift_vel", drift_vel_, + "The active drift velocity"); + drift_velocity_cmd.SetParameterName("drift_vel", true); + drift_velocity_cmd.SetUnitCategory("Velocity"); + + + G4GenericMessenger::Command& ELdrift_velocity_cmd = + msg_->DeclareProperty("ELdrift_vel", ELdrift_vel_, + "The EL drift velocity"); + ELdrift_velocity_cmd.SetParameterName("ELdrift_vel", true); + ELdrift_velocity_cmd.SetUnitCategory("Velocity"); + msg_->DeclareProperty("elfield", elfield_, "True if the EL field is on (full simulation), " "false if it's not (parametrized simulation."); @@ -329,7 +347,7 @@ void Next100FieldCage::BuildActive() G4double global_active_zpos = active_zpos_ - GetCoordOrigin().z(); field->SetCathodePosition(global_active_zpos + active_length_/2.); field->SetAnodePosition(global_active_zpos - active_length_/2.); - field->SetDriftVelocity(1. * mm/microsecond); + field->SetDriftVelocity(drift_vel_); field->SetTransverseDiffusion(drift_transv_diff_); field->SetLongitudinalDiffusion(drift_long_diff_); field->SetLifetime(e_lifetime_); @@ -721,10 +739,18 @@ void Next100FieldCage::BuildELRegion() G4double global_el_gap_zpos = el_gap_zpos_ - GetCoordOrigin().z(); el_field->SetCathodePosition(global_el_gap_zpos + el_gap_length_/2. + grid_thickn_); el_field->SetAnodePosition (global_el_gap_zpos - el_gap_length_/2. - grid_thickn_); - el_field->SetDriftVelocity(2.5 * mm/microsecond); + el_field->SetDriftVelocity(ELdrift_vel_); 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_)); + } + G4Region* el_region = new G4Region("EL_REGION"); el_region->SetUserInformation(el_field); el_region->AddRootLogicalVolume(el_gap_logic); diff --git a/source/geometries/Next100FieldCage.h b/source/geometries/Next100FieldCage.h index e964fe29a..135bf48a0 100644 --- a/source/geometries/Next100FieldCage.h +++ b/source/geometries/Next100FieldCage.h @@ -70,6 +70,10 @@ namespace nexus { G4double drift_transv_diff_, drift_long_diff_; G4double ELtransv_diff_; ///< transversal diffusion in the EL gap G4double ELlong_diff_; ///< longitudinal diffusion in the EL gap + + // Drift Velocities + G4double drift_vel_, ELdrift_vel_; + // Electric field G4bool elfield_; G4double ELelectric_field_; ///< electric field in the EL region diff --git a/source/geometries/Next100Vessel.cc b/source/geometries/Next100Vessel.cc index 463acf57b..783aa1f0f 100644 --- a/source/geometries/Next100Vessel.cc +++ b/source/geometries/Next100Vessel.cc @@ -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."); } + // 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_)); + } G4LogicalVolume* vessel_gas_logic = new G4LogicalVolume(vessel_gas_final_solid, vessel_gas_mat, "VESSEL_GAS"); diff --git a/source/materials/ArgonGasProperties.cc b/source/materials/ArgonGasProperties.cc index ed6c12c10..300a86b06 100644 --- a/source/materials/ArgonGasProperties.cc +++ b/source/materials/ArgonGasProperties.cc @@ -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; + + for (G4int i=0; i= 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; }