diff --git a/source/actions/LowMemoryTrackingAction.cc b/source/actions/LowMemoryTrackingAction.cc new file mode 100644 index 0000000000..4ab5bcd6a2 --- /dev/null +++ b/source/actions/LowMemoryTrackingAction.cc @@ -0,0 +1,98 @@ +// ---------------------------------------------------------------------------- +// nexus | LowMemoryTrackingAction.cc +// +// This class is a tracking action of the NEXT simulation. +// It stores in memory the trajectories of particles -- with the exception of +// optical photons, ionization electrons, and electrons starting in +// Rock or water -- with the relevant tracking information that will be +// saved to the output file. First version useful for TonneScale with no +// Cerenkov studies (June 2020). +// +// The NEXT Collaboration +// ---------------------------------------------------------------------------- + + +#include "LowMemoryTrackingAction.h" + +#include "Trajectory.h" +#include "TrajectoryMap.h" +#include "IonizationElectron.h" + +#include +#include +#include +#include +#include +#include + + + +using namespace nexus; + + + +LowMemoryTrackingAction::LowMemoryTrackingAction(): G4UserTrackingAction() +{ +} + + + +LowMemoryTrackingAction::~LowMemoryTrackingAction() +{ +} + + + +void LowMemoryTrackingAction::PreUserTrackingAction(const G4Track* track) +{ + // Do nothing if the track is an optical photon or an ionization electron + if (track->GetDefinition() == G4OpticalPhoton::Definition() || + track->GetDefinition() == IonizationElectron::Definition()) { + fpTrackingManager->SetStoreTrajectory(false); + return; + } + + if (track->GetParentID() == 0 || + track->GetParticleDefinition() == G4Neutron::Definition() || + track->GetParticleDefinition()->IsGeneralIon() ){ + + // Create a new trajectory associated to the track. + // N.B. If the processesing of a track is interrupted to be resumed + // later on (to process, for instance, its secondaries) more than + // one trajectory associated to the track will be created, but + // the event manager will merge them at some point. + G4VTrajectory* trj = new Trajectory(track); + + // Set the trajectory in the tracking manager + fpTrackingManager->SetStoreTrajectory(true); + fpTrackingManager->SetTrajectory(trj); + return; + } + + fpTrackingManager->SetStoreTrajectory(false); + } + + + +void LowMemoryTrackingAction::PostUserTrackingAction(const G4Track* track) +{ + // Do nothing if the track is an optical photon or an ionization electron + if (track->GetDefinition() == G4OpticalPhoton::Definition() || + track->GetDefinition() == IonizationElectron::Definition()) return; + + Trajectory* trj = (Trajectory*) TrajectoryMap::Get(track->GetTrackID()); + + // Do nothing if the track has no associated trajectory in the map + if (!trj) return; + + // Record final time and position of the track + trj->SetFinalPosition(track->GetPosition()); + trj->SetFinalTime(track->GetGlobalTime()); + trj->SetTrackLength(track->GetTrackLength()); + trj->SetFinalVolume(track->GetVolume()->GetName()); + trj->SetFinalMomentum(track->GetMomentum()); + + // Record last process of the track + G4String proc_name = track->GetStep()->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName(); + trj->SetFinalProcess(proc_name); +} diff --git a/source/actions/LowMemoryTrackingAction.h b/source/actions/LowMemoryTrackingAction.h new file mode 100644 index 0000000000..895ae3146c --- /dev/null +++ b/source/actions/LowMemoryTrackingAction.h @@ -0,0 +1,40 @@ +// ---------------------------------------------------------------------------- +// nexus | LowMemoryTrackingAction.h +// +// This class is a tracking action of the NEXT simulation. +// It stores in memory the trajectories of particles -- with the exception of +// optical photons, ionization electrons, and electrons starting in +// Rock or water -- with the relevant tracking information that will be +// saved to the output file. First version useful for TonneScale with no +// Cerenkov studies (June 2020). +// +// The NEXT Collaboration +// ---------------------------------------------------------------------------- + +#ifndef LOW_MEMORY_TRACKING_ACTION_H +#define LOW_MEMORY_TRACKING_ACTION_H + +#include + +class G4Track; + + +namespace nexus { + + // Low memory usage user tracking action + + class LowMemoryTrackingAction: public G4UserTrackingAction + { + public: + /// Constructor + LowMemoryTrackingAction(); + /// Destructor + virtual ~LowMemoryTrackingAction(); + + virtual void PreUserTrackingAction(const G4Track*); + virtual void PostUserTrackingAction(const G4Track*); + }; + +} + +#endif diff --git a/source/actions/NoSaveEventAction.cc b/source/actions/NoSaveEventAction.cc new file mode 100644 index 0000000000..47f75c494c --- /dev/null +++ b/source/actions/NoSaveEventAction.cc @@ -0,0 +1,35 @@ +#include "NoSaveEventAction.h" + +#include "PersistencyManager.h" + +#include + + +namespace nexus { + + + NoSaveEventAction::NoSaveEventAction(): + G4UserEventAction() + { + + } + + + + NoSaveEventAction::~NoSaveEventAction() + { + } + + + + void NoSaveEventAction::BeginOfEventAction(const G4Event* /*event*/) + { + PersistencyManager* pm = dynamic_cast + (G4VPersistencyManager::GetPersistencyManager()); + pm->StoreCurrentEvent(false); + } + + void NoSaveEventAction::EndOfEventAction(const G4Event* event) + { + } +} diff --git a/source/actions/NoSaveEventAction.h b/source/actions/NoSaveEventAction.h new file mode 100644 index 0000000000..ca95cedca9 --- /dev/null +++ b/source/actions/NoSaveEventAction.h @@ -0,0 +1,39 @@ +// ---------------------------------------------------------------------------- +// nexus | NoSaveEventAction.h +// +// This event action saves no event information to file +// +// The NEXT Collaboration +// ---------------------------------------------------------------------------- + +#ifndef NOSAVE_EVENT_ACTION_H +#define NOSAVE_EVENT_ACTION_H + +#include +#include + +class G4Event; + +namespace nexus { + + /// This class is an Event Action which forces nexus + /// Not to store any event information + + class NoSaveEventAction: public G4UserEventAction + { + public: + /// Constructor + NoSaveEventAction(); + /// Destructor + ~NoSaveEventAction(); + + /// Hook at the beginning of the event loop + void BeginOfEventAction(const G4Event*); + /// Hook at the end of the event loop + void EndOfEventAction(const G4Event*); + + }; + +} // namespace nexus + +#endif diff --git a/source/base/ActionsFactory.cc b/source/base/ActionsFactory.cc index 0e31107614..e9acc7f0fd 100644 --- a/source/base/ActionsFactory.cc +++ b/source/base/ActionsFactory.cc @@ -58,6 +58,7 @@ G4UserRunAction* ActionsFactory::CreateRunAction() const #include "DefaultEventAction.h" #include "SaveAllEventAction.h" #include "MuonsEventAction.h" +#include "NoSaveEventAction.h" G4UserEventAction* ActionsFactory::CreateEventAction() const { @@ -69,6 +70,7 @@ G4UserEventAction* ActionsFactory::CreateEventAction() const else if (evtact_name_ == "MUONS") p = new MuonsEventAction(); + else if (evtact_name_ == "NO_SAVE") p = new NoSaveEventAction(); else { G4String err = "Unknown user event action: " + evtact_name_; G4Exception("[ActionsFactory]", "CreateEventAction()", JustWarning, err); @@ -82,6 +84,7 @@ G4UserEventAction* ActionsFactory::CreateEventAction() const #include "DefaultTrackingAction.h" #include "ValidationTrackingAction.h" #include "OpticalTrackingAction.h" +#include "LowMemoryTrackingAction.h" G4UserTrackingAction* ActionsFactory::CreateTrackingAction() const { @@ -93,6 +96,8 @@ G4UserTrackingAction* ActionsFactory::CreateTrackingAction() const else if (trkact_name_ == "OPTICAL") p = new OpticalTrackingAction(); + else if (trkact_name_ == "LOW_MEMORY") p = new LowMemoryTrackingAction(); + else { G4String err = "Unknown user tracking action: " + trkact_name_; G4Exception("[ActionsFactory]", "CreateTrackingAction()", diff --git a/source/geometries/NextTonScale.cc b/source/geometries/NextTonScale.cc index 79b062c058..02bb5a006f 100644 --- a/source/geometries/NextTonScale.cc +++ b/source/geometries/NextTonScale.cc @@ -24,6 +24,8 @@ #include #include #include +#include +#include using namespace nexus; @@ -33,6 +35,7 @@ NextTonScale::NextTonScale(): msg_(nullptr), gas_("naturalXe"), gas_pressure_(15.*bar), gas_temperature_(300.*kelvin), detector_diam_(0.), detector_length_(0.), + rock_thickn_(2.*m), rock_diam_(0.), tank_size_(0.), tank_thickn_(1.*cm), water_thickn_(3.*m), vessel_thickn_(1.5*mm), ics_thickn_(12.*cm), endcap_hollow_(20.*cm), fcage_thickn_(1.*cm), @@ -73,12 +76,13 @@ void NextTonScale::Construct() 2.*ics_thickn_ + 2.*endcap_hollow_ + 2.*vessel_thickn_; tank_size_ = std::max(detector_diam_, detector_length_) + 2.*water_thickn_ + 2.*tank_thickn_; + rock_diam_ = tank_size_ + 2. * cm + 2. *m; // LABORATORY //////////////////////////////////////////// // Volume of air that contains all other detector volumes. G4String lab_name = "LABORATORY"; - G4double lab_size = tank_size_ + 4.*m; // Clearance of 2 m around water tank + G4double lab_size = rock_diam_ + 6.*m; // Clearance of 2 m around water tank and rock G4Material* lab_material = G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR"); @@ -96,12 +100,57 @@ void NextTonScale::Construct() G4LogicalVolume* mother_logic_vol = lab_logic_vol; + mother_logic_vol = ConstructLabRock(mother_logic_vol); mother_logic_vol = ConstructWaterTank(mother_logic_vol); mother_logic_vol = ConstructVesselAndICS(mother_logic_vol); mother_logic_vol = ConstructFieldCageAndReadout(mother_logic_vol); } +G4LogicalVolume* NextTonScale::ConstructLabRock(G4LogicalVolume* mother_logic_vol) +{ + // Lab rock // + // Tube of rock (Standard or G4_CONCRETE) with additional cylinder above tank + G4String top_rock_name = "LAB_ROCK_TOP"; + G4String side_rock_name = "LAB_ROCK_SURROUND"; + G4Material* rock_material = MaterialsList::StandardRock(); + + G4Tubs* top_rock_solid = + new G4Tubs(top_rock_name, 0, rock_diam_ / 2., + rock_thickn_ / 2., 0., 360. * deg); + G4double side_rock_height = tank_size_ + 1. * cm; + G4double side_rock_inner_diam = tank_size_ + 2. * cm; + G4Tubs* side_rock_solid = + new G4Tubs(side_rock_name, side_rock_inner_diam / 2., rock_diam_ / 2., + side_rock_height / 2., 0., 360. * deg); + + G4Region *rock_region = new G4Region("LAB_ROCK"); + G4ProductionCuts *rock_prod_cut = new G4ProductionCuts(); + rock_prod_cut->SetProductionCut(50. * cm); + G4LogicalVolume* top_rock_logic = + new G4LogicalVolume(top_rock_solid, rock_material, top_rock_name); + rock_region->AddRootLogicalVolume(top_rock_logic); + + G4LogicalVolume* side_rock_logic = + new G4LogicalVolume(side_rock_solid, rock_material, side_rock_name); + rock_region->AddRootLogicalVolume(side_rock_logic); + rock_region->SetProductionCuts(rock_prod_cut); + + G4RotationMatrix* rock_rotation = new G4RotationMatrix(); + rock_rotation->rotateX(90. * deg); + G4double top_rock_y = tank_size_ / 2. + 1.*cm + rock_thickn_ / 2; + G4ThreeVector top_rock_pos(0., top_rock_y, 0.); + new G4PVPlacement(rock_rotation, top_rock_pos, top_rock_logic, + top_rock_name, mother_logic_vol, false, 0, true); + + G4double side_rock_y = 0.5 * cm; + new G4PVPlacement(rock_rotation, G4ThreeVector(0., side_rock_y, 0.), + side_rock_logic, side_rock_name, mother_logic_vol, + false, 0, true); + return mother_logic_vol; +} + + G4LogicalVolume* NextTonScale::ConstructWaterTank(G4LogicalVolume* mother_logic_vol) { // WATER TANK //////////////////////////////////////////// @@ -146,9 +195,10 @@ G4LogicalVolume* NextTonScale::ConstructWaterTank(G4LogicalVolume* mother_logic_ ////////////////////////////////////////////////////////// - muon_gen_ = new MuonsPointSampler(tank_size_/2. + 50. * cm, - tank_size_/2. + 1. * cm, - tank_size_/2. + 50. * cm); + G4double mu_gen_x = rock_diam_ / 2.; + G4double mu_gen_y = tank_size_/2. + rock_thickn_ + 2. * cm; + G4double mu_gen_z = rock_diam_ / 2.; + muon_gen_ = new MuonsPointSampler(mu_gen_x, mu_gen_y, mu_gen_z, true); // Primarily for neutron generation external_gen_ = new CylinderPointSampler(tank_size_/2. + 1 * cm, @@ -482,6 +532,13 @@ void NextTonScale::DefineConfigurationParameters() gas_temperature_cmd.SetParameterName("gas_temperature", false); gas_temperature_cmd.SetRange("gas_temperature>0."); + G4GenericMessenger::Command& rock_thickn_cmd = + msg_->DeclareProperty("rock_thickness", rock_thickn_, + "Thickness of surrounding rock to simulate."); + rock_thickn_cmd.SetUnitCategory("Length"); + rock_thickn_cmd.SetParameterName("rock_thickness", false); + rock_thickn_cmd.SetRange("rock_thickness>=0."); + G4GenericMessenger::Command& active_diam_cmd = msg_->DeclareProperty("active_diam", active_diam_, "Diameter of the active diameter of the detector."); diff --git a/source/geometries/NextTonScale.h b/source/geometries/NextTonScale.h index feb02de7dc..e56ede13db 100644 --- a/source/geometries/NextTonScale.h +++ b/source/geometries/NextTonScale.h @@ -34,6 +34,7 @@ namespace nexus { private: // + G4LogicalVolume* ConstructLabRock(G4LogicalVolume*); G4LogicalVolume* ConstructWaterTank(G4LogicalVolume*); G4LogicalVolume* ConstructVesselAndICS(G4LogicalVolume*); G4LogicalVolume* ConstructFieldCageAndReadout(G4LogicalVolume*); @@ -50,6 +51,7 @@ namespace nexus { G4String gas_; G4double gas_pressure_, gas_temperature_; G4double detector_diam_, detector_length_; + G4double rock_thickn_, rock_diam_; G4double tank_size_, tank_thickn_, water_thickn_; G4double vessel_thickn_, ics_thickn_, endcap_hollow_; G4double fcage_thickn_, active_diam_, active_length_; diff --git a/source/materials/MaterialsList.cc b/source/materials/MaterialsList.cc index 7a31d6205f..e483900ab4 100644 --- a/source/materials/MaterialsList.cc +++ b/source/materials/MaterialsList.cc @@ -265,6 +265,21 @@ G4Material* MaterialsList::GXeHe(G4double pressure, } +G4Material* MaterialsList::StandardRock(){ + + G4String name = "StandardRock"; + + G4Material* mat = G4Material::GetMaterial(name, false); + + if (mat == 0){ + G4double Z = 11; + G4double A = 22; + G4double rho = 2.65 * g / cm3; + mat = new G4Material(name, Z, A, rho, kStateSolid); + } + return mat; +} + diff --git a/source/materials/MaterialsList.h b/source/materials/MaterialsList.h index 83557c1583..ab38eef3d7 100644 --- a/source/materials/MaterialsList.h +++ b/source/materials/MaterialsList.h @@ -134,6 +134,9 @@ namespace nexus { static G4Material* CopyMaterial(G4Material*, const G4String&); + // Standard rock + static G4Material* StandardRock(); + private: /// Constructor (hidden) MaterialsList(); diff --git a/source/utils/MuonsPointSampler.cc b/source/utils/MuonsPointSampler.cc index a34cc8e53d..975e440363 100644 --- a/source/utils/MuonsPointSampler.cc +++ b/source/utils/MuonsPointSampler.cc @@ -21,8 +21,8 @@ namespace nexus { using namespace CLHEP; MuonsPointSampler::MuonsPointSampler - (G4double x, G4double yPoint, G4double z): - x_(x),yPoint_(yPoint),z_(z) + (G4double x, G4double yPoint, G4double z, G4bool disc): + x_(x),yPoint_(yPoint),z_(z),disc_(disc) { } @@ -42,6 +42,12 @@ namespace nexus { // y is fixed G4double x = -x_ + G4UniformRand()*2*x_; G4double z = -z_ + G4UniformRand()*2*z_; + if (disc_){ + while (x*x / (x_*x_) + z*z / (z_*z_) > 1){ + x = -x_ + G4UniformRand()*2*x_; + z = -z_ + G4UniformRand()*2*z_; + } + } G4ThreeVector mypoint(x, yPoint_, z); diff --git a/source/utils/MuonsPointSampler.h b/source/utils/MuonsPointSampler.h index a2dc8fd129..f93bc68829 100644 --- a/source/utils/MuonsPointSampler.h +++ b/source/utils/MuonsPointSampler.h @@ -20,8 +20,7 @@ namespace nexus { public: /// Constructor - MuonsPointSampler(G4double x, G4double yPoint, G4double z); - + MuonsPointSampler(G4double x, G4double yPoint, G4double z, G4bool disc=false); /// Destructor ~MuonsPointSampler(); @@ -37,6 +36,7 @@ namespace nexus { G4ThreeVector GetXZPointInMuonsPlane(); G4double x_, yPoint_,z_; + G4bool disc_; };