From 1eb47e375bfb119a722b8bfb0209c6b798a90ae5 Mon Sep 17 00:00:00 2001 From: Andrew Laing Date: Fri, 20 Nov 2020 20:28:57 +0100 Subject: [PATCH 01/25] Add ProjectToRegion to BaseGeometry --- source/geometries/GeometryBase.h | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/source/geometries/GeometryBase.h b/source/geometries/GeometryBase.h index 94665e4b3e..6a70b4c952 100644 --- a/source/geometries/GeometryBase.h +++ b/source/geometries/GeometryBase.h @@ -34,6 +34,12 @@ namespace nexus { /// Returns a point within a given region of the geometry virtual G4ThreeVector GenerateVertex(const G4String&) const; + /// Returns a point within a region projecting from a + /// given point backwards along a line. + virtual G4ThreeVector ProjectToRegion(const G4String&, + const G4ThreeVector&, + const G4ThreeVector&) const; + /// Returns the span (maximum dimension) of the geometry G4double GetSpan(); @@ -102,7 +108,12 @@ namespace nexus { inline G4ThreeVector GeometryBase::GenerateVertex(const G4String&) const { return G4ThreeVector(0., 0., 0.); } - inline void GeometryBase::SetSpan(G4double s) { span_ = s; } + inline G4ThreeVector BaseGeometry::ProjectToRegion(const G4String&, + const G4ThreeVector&, + const G4ThreeVector&) const + { return G4ThreeVector(0., 0., 0.); } + + inline void BaseGeometry::SetSpan(G4double s) { span_ = s; } inline G4double GeometryBase::GetSpan() { return span_; } From 9f623b0f442cef32cff295594b620bb57258d74a Mon Sep 17 00:00:00 2001 From: Andrew Laing Date: Fri, 20 Nov 2020 20:29:36 +0100 Subject: [PATCH 02/25] Add function to find the intersection of a ray with BoxPointSampler --- source/utils/BoxPointSampler.cc | 19 +++++++++++++++++++ source/utils/BoxPointSampler.h | 4 ++++ 2 files changed, 23 insertions(+) diff --git a/source/utils/BoxPointSampler.cc b/source/utils/BoxPointSampler.cc index 5e3cf2e5e3..dc300bfa49 100644 --- a/source/utils/BoxPointSampler.cc +++ b/source/utils/BoxPointSampler.cc @@ -206,4 +206,23 @@ namespace nexus { return real_pos; } + G4ThreeVector BoxPointSampler::GetIntersect(const G4ThreeVector& point, + const G4ThreeVector& dir) + { + // Get the +ve movements to intersect with the faces of the box. + G4double tx = (-inner_x_ - point.x()) / dir.x(); + if (tx < 0) tx = (inner_x_ - point.x()) / dir.x(); + + G4double ty = (-inner_y_ - point.y()) / dir.y(); + if (ty < 0) ty = (inner_y_ - point.y()) / dir.y(); + + G4double tz = (-inner_z_ - point.z()) / dir.z(); + if (tz < 0) tz = (inner_z_ - point.z()) / dir.z(); + + // The minimum of the tx, ty, tz gives the intersection point. + G4double tmin = std::min({tx, ty, tz}); + + return RotateAndTranslate(point + tmin * dir); + } + } // end namespace nexus diff --git a/source/utils/BoxPointSampler.h b/source/utils/BoxPointSampler.h index fb44196f93..84e2349c52 100644 --- a/source/utils/BoxPointSampler.h +++ b/source/utils/BoxPointSampler.h @@ -32,6 +32,10 @@ namespace nexus { /// Return vertex within region of the chamber G4ThreeVector GenerateVertex(const G4String& region); + /// Return the intersect point along dir + G4ThreeVector GetIntersect(const G4ThreeVector& point, + const G4ThreeVector& dir); + private: G4double GetLength(G4double origin, G4double max_length); G4ThreeVector RotateAndTranslate(G4ThreeVector position); From d57632010ec3a82b8b1feaee021cfbb2432b18d1 Mon Sep 17 00:00:00 2001 From: Andrew Laing Date: Fri, 20 Nov 2020 20:30:10 +0100 Subject: [PATCH 03/25] Add Ray intersection calculation to CylinderPointSampler2020 --- source/utils/CylinderPointSampler2020.cc | 40 ++++++++++++++++++++++++ source/utils/CylinderPointSampler2020.h | 4 +++ 2 files changed, 44 insertions(+) diff --git a/source/utils/CylinderPointSampler2020.cc b/source/utils/CylinderPointSampler2020.cc index 83bd1fbf3a..01360031ce 100644 --- a/source/utils/CylinderPointSampler2020.cc +++ b/source/utils/CylinderPointSampler2020.cc @@ -104,6 +104,46 @@ namespace nexus { } + G4ThreeVector + CylinderPointSampler2020::GetIntersect(const G4ThreeVector& point, + const G4ThreeVector& dir) + { + // Projects for +ve direction (for backwards send -ve dir) + // To find the intersection with the Cylinder. + // First we want to solve for the ray intersection with + // an infinite barrel, so where point + tdir satisfies + // x^2 + y^2 - r^2 = 0 (valid for all z before rotation). + // Solve the quadratic equation and take the +ve t + G4double a = dir.x()*dir.x() + dir.y()*dir.y(); + G4double b = 2 * (point.x()*dir.x() + point.y()*dir.y()); + G4double c = point.x()*point.x() + point.y()*point.y() - minRad_*minRad_; + G4double determinant = b*b - 4*a*c; + if (determinant < 0) + // No intersection return origin. + return RotateAndTranslate(G4ThreeVector(0., 0., 0.)); + + G4double t = -b + std::sqrt(determinant) / (2 * a); + if (t < 0) t = -b - std::sqrt(determinant) / (2 * a); + // If still -ve we're outside the volume giv origin + if (t < 0) + return RotateAndTranslate(G4ThreeVector(0., 0., 0.)); + + G4ThreeVector barrel_intersect = point + t * dir; + + // If the z of the prediction is within the barrel limits we're done. + if (std::abs(barrel_intersect.z()) <= halfLength_) + return RotateAndTranslate(barrel_intersect); + + // If not, we need to check the endcaps. + if (barrel_intersect.z() < 0) + // We're beyond the -ve endcap. + t = (-halfLength_ - point.z()) / dir.z(); + else + t = (halfLength_ - point.z()) / dir.z(); + return RotateAndTranslate(point + t * dir); + } + + G4double CylinderPointSampler2020::GetRadius(G4double innerRad, G4double outerRad) { diff --git a/source/utils/CylinderPointSampler2020.h b/source/utils/CylinderPointSampler2020.h index af2d0ff5b8..d7e116aac5 100644 --- a/source/utils/CylinderPointSampler2020.h +++ b/source/utils/CylinderPointSampler2020.h @@ -41,6 +41,10 @@ namespace nexus { // Returns vertex within region of the chamber G4ThreeVector GenerateVertex(const G4String& region); + /// Return the intersect point along dir + G4ThreeVector GetIntersect(const G4ThreeVector& point, + const G4ThreeVector& dir); + private: G4double GetRadius(G4double innerRad, G4double outerRad); G4double GetPhi(); From d6310f309b63edabcfc525aeaf5f19f2599fdf2f Mon Sep 17 00:00:00 2001 From: Paola Ferrario Date: Thu, 28 Sep 2023 09:31:24 +0200 Subject: [PATCH 04/25] Add machinery to use Projections for vertex generation Adds to NEW and Next100 for the EXTERNAL generator of the lead shielding --- source/geometries/Next100.cc | 19 +++++++++++++++++++ source/geometries/Next100.h | 6 ++++++ source/geometries/Next100Shielding.cc | 27 +++++++++++++++++++++++++++ source/geometries/Next100Shielding.h | 6 ++++++ source/geometries/NextNew.cc | 18 ++++++++++++++++++ source/geometries/NextNew.h | 6 ++++++ 6 files changed, 82 insertions(+) diff --git a/source/geometries/Next100.cc b/source/geometries/Next100.cc index 80c7471882..446c5ed1dd 100644 --- a/source/geometries/Next100.cc +++ b/source/geometries/Next100.cc @@ -271,4 +271,23 @@ namespace nexus { return vertex; } + + G4ThreeVector Next100::ProjectToRegion(const G4String& region, + const G4ThreeVector& point, + const G4ThreeVector& dir) const + { + // Project backwards along dir from point to find the first intersection + // with region. + G4ThreeVector vertex(0., 0., 0.); + if (region == "EXTERNAL"){ + return shielding_->ProjectToRegion(region, point, dir); + } + else { + G4Exception("[Next100]", "ProjectToRegion()", FatalException, + "Unknown vertex generation region!"); + } + + return vertex + G4ThreeVector(0., 0., -gate_zpos_in_vessel_); + } + } //end namespace nexus diff --git a/source/geometries/Next100.h b/source/geometries/Next100.h index 828668d619..11496d64db 100644 --- a/source/geometries/Next100.h +++ b/source/geometries/Next100.h @@ -38,6 +38,12 @@ namespace nexus { /// Generate a vertex within a given region of the geometry G4ThreeVector GenerateVertex(const G4String& region) const; + /// Returns a point within a region projecting from a + /// given point backwards along a line. + G4ThreeVector ProjectToRegion(const G4String& region, + const G4ThreeVector& point, + const G4ThreeVector& dir) const; + private: void BuildLab(); diff --git a/source/geometries/Next100Shielding.cc b/source/geometries/Next100Shielding.cc index 7a12a6932d..18c87de068 100644 --- a/source/geometries/Next100Shielding.cc +++ b/source/geometries/Next100Shielding.cc @@ -371,6 +371,13 @@ namespace nexus { lead_gen_ = new BoxPointSampler(steel_x, steel_y, steel_z, 5.*cm, G4ThreeVector(0., 0., 0.), 0); + G4double shield_diag = std::sqrt(lead_x_*lead_x_ + lead_y_*lead_y_ + lead_z_*lead_z_); + G4double ext_offset = 1. * cm; + external_gen_ = new BoxPointSampler(shield_diag / 2. + ext_offset, + shield_diag / 2. + ext_offset, + shield_diag / 2. + ext_offset, + 1. * mm, G4ThreeVector(0.,0.,0.), 0); + steel_gen_ = new BoxPointSampler(shield_x_, shield_y_, shield_z_, steel_thickness_, G4ThreeVector(0., -beam_thickness_2/2., 0.), 0); @@ -787,4 +794,24 @@ namespace nexus { return vertex; } + + + G4ThreeVector Next100Shielding::ProjectToRegion(const G4String& region, + const G4ThreeVector& point, + const G4ThreeVector& dir) const + { + // Project backwards along dir from point to find the first intersection + // with region. + G4ThreeVector vertex(0., 0., 0.); + if (region == "EXTERNAL"){ + return external_gen_->GetIntersect(point, dir); + } + else { + G4Exception("[Next100Shielding]", "ProjectToRegion()", FatalException, + "Unknown vertex generation region!"); + } + + return vertex; + } + } //end namespace nexus diff --git a/source/geometries/Next100Shielding.h b/source/geometries/Next100Shielding.h index 1110865599..db1996b9d1 100644 --- a/source/geometries/Next100Shielding.h +++ b/source/geometries/Next100Shielding.h @@ -41,6 +41,12 @@ namespace nexus { G4double GetHeight() const; + /// Returns a point within a region projecting from a + /// given point backwards along a line. + G4ThreeVector ProjectToRegion(const G4String& region, + const G4ThreeVector& point, + const G4ThreeVector& dir) const; + /// Builder void Construct(); diff --git a/source/geometries/NextNew.cc b/source/geometries/NextNew.cc index 261f3537c2..e33e3df730 100644 --- a/source/geometries/NextNew.cc +++ b/source/geometries/NextNew.cc @@ -599,5 +599,23 @@ namespace nexus { return vertex; } + G4ThreeVector NextNew::ProjectToRegion(const G4String& region, + const G4ThreeVector& point, + const G4ThreeVector& dir) const + { + // Project backwards along dir from point to find the first intersection + // with region. + G4ThreeVector vertex(0., 0., 0.); + if (region == "EXTERNAL"){ + return shielding_->ProjectToRegion(region, point, dir); + } + else { + G4Exception("[NextNew]", "ProjectToRegion()", FatalException, + "Unknown vertex generation region!"); + } + + return vertex + displ_; + } + } //end namespace nexus diff --git a/source/geometries/NextNew.h b/source/geometries/NextNew.h index 45858903ef..d1029458d7 100644 --- a/source/geometries/NextNew.h +++ b/source/geometries/NextNew.h @@ -46,6 +46,12 @@ namespace nexus { /// Generate a vertex within a given region of the geometry G4ThreeVector GenerateVertex(const G4String& region) const; + /// Returns a point within a region projecting from a + /// given point backwards along a line. + G4ThreeVector ProjectToRegion(const G4String& region, + const G4ThreeVector& point, + const G4ThreeVector& dir) const; + private: void BuildExtScintillator(G4ThreeVector pos, const G4RotationMatrix& rot); void Construct(); From a1595b022f87708887a3110c3a1686d1a7521b40 Mon Sep 17 00:00:00 2001 From: Paola Ferrario Date: Thu, 28 Sep 2023 09:46:03 +0200 Subject: [PATCH 05/25] Update MuonAngleGenerator to use a rotatic disc and projection Points generated on disc then ray projected to the requested region. --- source/generators/MuonAngleGenerator.cc | 202 ++++++++++++++++++++++++ source/generators/MuonAngleGenerator.h | 82 ++++++++++ 2 files changed, 284 insertions(+) create mode 100644 source/generators/MuonAngleGenerator.cc create mode 100644 source/generators/MuonAngleGenerator.h diff --git a/source/generators/MuonAngleGenerator.cc b/source/generators/MuonAngleGenerator.cc new file mode 100644 index 0000000000..6dacde506b --- /dev/null +++ b/source/generators/MuonAngleGenerator.cc @@ -0,0 +1,202 @@ +// ---------------------------------------------------------------------------- +// nexus | MuonAngleGenerator.cc +// +// This class is the primary generator of muons following an angular +// distribution measured in the LSC. +// +// The NEXT Collaboration +// ---------------------------------------------------------------------------- + +#include "MuonAngleGenerator.h" +#include "DetectorConstruction.h" +#include "BaseGeometry.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "MuonsPointSampler.h" +#include "AddUserInfoToPV.h" + +#include +#include "TFile.h" +#include "TH2F.h" +#include "CLHEP/Units/SystemOfUnits.h" + +using namespace nexus; +using namespace CLHEP; + + +MuonAngleGenerator::MuonAngleGenerator(): + G4VPrimaryGenerator(), msg_(0), particle_definition_(0), + angular_generation_(true), rPhi_(NULL), energy_min_(0.), + energy_max_(0.), gen_rad_(223.33*cm), distribution_(0), geom_(0)//, geom_solid_(0) +{ + msg_ = new G4GenericMessenger(this, "/Generator/MuonAngleGenerator/", + "Control commands of muongenerator."); + + G4GenericMessenger::Command& min_energy = + msg_->DeclareProperty("min_energy", energy_min_, "Set minimum kinetic energy of the particle."); + min_energy.SetUnitCategory("Energy"); + min_energy.SetParameterName("min_energy", false); + min_energy.SetRange("min_energy>0."); + + G4GenericMessenger::Command& max_energy = + msg_->DeclareProperty("max_energy", energy_max_, "Set maximum kinetic energy of the particle"); + max_energy.SetUnitCategory("Energy"); + max_energy.SetParameterName("max_energy", false); + max_energy.SetRange("max_energy>0."); + + // defaults to length of the NEXT100 shielding diagonal if not present. + G4GenericMessenger::Command& generation_radius = + msg_->DeclareProperty("gen_rad", gen_rad_, "Set radius for generation disc"); + generation_radius.SetUnitCategory("Length"); + generation_radius.SetParameterName("gen_rad", false); + generation_radius.SetRange("gen_rad>0."); + + msg_->DeclareProperty("region", region_, + "Set the region of the geometry where the vertex will be generated."); + + msg_->DeclareProperty("angles_on", angular_generation_, + "Distribute muon directions according to file?"); + + msg_->DeclareProperty("angle_file", ang_file_, + "Name of the file containing angular distribution."); + msg_->DeclareProperty("angle_dist", dist_name_, + "Name of the angular distribution histogram."); + + G4GenericMessenger::Command& rotation = + msg_->DeclareProperty("azimuth_rotation", axis_rotation_, + "Angle between north and nexus z in anticlockwise"); + rotation.SetUnitCategory("Angle"); + rotation.SetParameterName("azimuth", false); + rotation.SetRange("azimuth>0."); + + DetectorConstruction* detconst = (DetectorConstruction*) G4RunManager::GetRunManager()->GetUserDetectorConstruction(); + geom_ = detconst->GetGeometry(); + +} + + +MuonAngleGenerator::~MuonAngleGenerator() +{ + + delete msg_; +} + + +void MuonAngleGenerator::SetupAngles() +{ + // Rotation from the axes used in file. + // Rotates anticlockwise about Y. + rPhi_ = new G4RotationMatrix(); + rPhi_->rotateY(-axis_rotation_); + + // Get the Angular distribution from file. + TFile angle_file(ang_file_); + angle_file.GetObject(dist_name_, distribution_); + distribution_->SetDirectory(0); + angle_file.Close(); + +} + + +void MuonAngleGenerator::GeneratePrimaryVertex(G4Event* event) +{ + + if (angular_generation_ && rPhi_ == NULL) + SetupAngles(); + + particle_definition_ = + G4ParticleTable::GetParticleTable()->FindParticle(MuonCharge()); + if (!particle_definition_) + G4Exception("[MuonAngleGenerator]", "SetParticleDefinition()", + FatalException, " can not create a muon "); + + // Generate uniform random energy in [E_min, E_max] + G4double kinetic_energy = RandomEnergy(); + + // Particle propierties + G4double mass = particle_definition_->GetPDGMass(); + G4double energy = kinetic_energy + mass; + G4double pmod = std::sqrt(energy*energy - mass*mass); + + G4ThreeVector p_dir(0., -1., 0.); + if (angular_generation_) + GetDirection(p_dir); + G4ThreeVector position = ProjectToVertex(p_dir); + G4double px = pmod * p_dir.x(); + G4double py = pmod * p_dir.y(); + G4double pz = pmod * p_dir.z(); + + // Particle generated at start-of-event + G4double time = 0.; + // Create a new vertex + G4PrimaryVertex* vertex = new G4PrimaryVertex(position, time); + + // Create the new primary particle and set it some properties + G4PrimaryParticle* particle = + new G4PrimaryParticle(particle_definition_, px, py, pz); + + // Add particle to the vertex and this to the event + vertex->SetPrimary(particle); + event->AddPrimaryVertex(vertex); +} + + +G4double MuonAngleGenerator::RandomEnergy() const +{ + if (energy_max_ == energy_min_) + return energy_min_; + else + return (G4UniformRand()*(energy_max_ - energy_min_) + energy_min_); +} + + +G4String MuonAngleGenerator::MuonCharge() const +{ + G4double rndCh = 2.3 *G4UniformRand(); //From PDG cosmic muons mu+/mu- = 1.3 + if (rndCh <1.3) + return "mu+"; + else + return "mu-"; +} + + +void MuonAngleGenerator::GetDirection(G4ThreeVector& dir) +{ + // GetAngles from file?? Azimuth defined anticlockwise + // From north + G4double zenith = 0.; + G4double azimuth = 0.; + distribution_->GetRandom2(azimuth, zenith); + // !! Current distribution in units of pi + zenith *= pi; + azimuth *= pi; + + dir.setX(sin(zenith) * sin(azimuth)); + dir.setY(-cos(zenith)); + dir.setZ(-sin(zenith) * cos(azimuth)); + + dir *= *rPhi_; +} + + +G4ThreeVector MuonAngleGenerator::ProjectToVertex(const G4ThreeVector& dir) +{ + // Postion in disc (need to sort size, member function) + G4double rad = gen_rad_ * std::sqrt(G4UniformRand()); + G4double ang = 2 * G4UniformRand() * pi; + + // Continue assuming that Y is vertical and z drift, + // valid for NEW and NEXT-100 (at least) + G4ThreeVector point(rad * std::cos(ang), 0., rad * std::sin(ang)); + point.rotate(pi / 2 - dir.angle(point), dir.cross(point)); + + // Now project back to the requested region intersection. + return geom_->ProjectToRegion(region_, point, -dir); +} diff --git a/source/generators/MuonAngleGenerator.h b/source/generators/MuonAngleGenerator.h new file mode 100644 index 0000000000..cb48f229a0 --- /dev/null +++ b/source/generators/MuonAngleGenerator.h @@ -0,0 +1,82 @@ +// ---------------------------------------------------------------------------- +// nexus | MuonAngleGenerator.h +// +// This class is the primary generator of muons following an angular +// distribution measured in the LSC. +// +// The NEXT Collaboration +// ---------------------------------------------------------------------------- + +#ifndef MUON_ANGLE_GENERATOR_H +#define MUON_ANGLE_GENERATOR_H + +#include +#include + +class G4GenericMessenger; +class G4Event; +class G4ParticleDefinition; +//class G4VSolid; + +class TH2F; + + +namespace nexus { + + class BaseGeometry; + + class MuonAngleGenerator: public G4VPrimaryGenerator + { + public: + /// Constructor + MuonAngleGenerator(); + /// Destructor + ~MuonAngleGenerator(); + + /// This method is invoked at the beginning of the event. It sets + /// a primary vertex (that is, a particle in a given position and time) + /// in the event. + void GeneratePrimaryVertex(G4Event*); + + private: + + // Sets the rotation angle and the spectra to + // be read for angle generation as well as + // setting the overlap volume for filtering. + void SetupAngles(); + + /// Generate a random kinetic energy with flat probability in + // the interval [energy_min, energy_max]. + G4double RandomEnergy() const; + G4String MuonCharge() const; + + void GetDirection(G4ThreeVector& dir); + + G4ThreeVector ProjectToVertex(const G4ThreeVector& dir); + + private: + G4GenericMessenger* msg_; + + G4ParticleDefinition* particle_definition_; + + G4bool angular_generation_; ///< Distribution or all downwards + G4double axis_rotation_; ///< Angle between North and +z + G4RotationMatrix *rPhi_; ///< Rotation to adjust axes + + G4double energy_min_; ///< Minimum kinetic energy + G4double energy_max_; ///< Maximum kinetic energy + + G4String region_; ///< Name of generator region + G4String ang_file_; ///< Name of file with distributions + G4String dist_name_; ///< Name of distribution in file + G4double gen_rad_; ///< Radius of disc for generation + + TH2F * distribution_; ///< Anglular distribution + + const BaseGeometry* geom_; ///< Pointer to the detector geometry + + }; + +} // end namespace nexus + +#endif From fceb48654f05c05831fa81e33d4cec737d46bd71 Mon Sep 17 00:00:00 2001 From: Andrew Laing Date: Tue, 24 Nov 2020 18:31:58 +0100 Subject: [PATCH 06/25] Take into account arbitrary rotation and translation in BoxSampler --- source/utils/BoxPointSampler.cc | 35 ++++++++++++++++++++++++++------- source/utils/BoxPointSampler.h | 1 + 2 files changed, 29 insertions(+), 7 deletions(-) diff --git a/source/utils/BoxPointSampler.cc b/source/utils/BoxPointSampler.cc index dc300bfa49..eac5d2f5c9 100644 --- a/source/utils/BoxPointSampler.cc +++ b/source/utils/BoxPointSampler.cc @@ -206,23 +206,44 @@ namespace nexus { return real_pos; } + void BoxPointSampler::InvertRotationAndTranslation(G4ThreeVector& vec, + bool translate) + { + if (translate) + vec -= origin_; + if (rotation_) + vec.rotate(-rotation_->delta(), rotation_->axis()); + } + G4ThreeVector BoxPointSampler::GetIntersect(const G4ThreeVector& point, const G4ThreeVector& dir) { + // The point and direction should be in the lab frame. + // Need to rotate into the Sampler frame to get the intersection + // then rotate back into lab frame. + G4ThreeVector local_point = point; + InvertRotationAndTranslation(local_point); + G4ThreeVector local_dir = dir; + InvertRotationAndTranslation(local_dir, false); + // Get the +ve movements to intersect with the faces of the box. - G4double tx = (-inner_x_ - point.x()) / dir.x(); - if (tx < 0) tx = (inner_x_ - point.x()) / dir.x(); + G4double tx = (-inner_x_ - local_point.x()) / local_dir.x(); + if (tx < 0) tx = (inner_x_ - local_point.x()) / local_dir.x(); - G4double ty = (-inner_y_ - point.y()) / dir.y(); - if (ty < 0) ty = (inner_y_ - point.y()) / dir.y(); + G4double ty = (-inner_y_ - local_point.y()) / local_dir.y(); + if (ty < 0) ty = (inner_y_ - local_point.y()) / local_dir.y(); - G4double tz = (-inner_z_ - point.z()) / dir.z(); - if (tz < 0) tz = (inner_z_ - point.z()) / dir.z(); + G4double tz = (-inner_z_ - local_point.z()) / local_dir.z(); + if (tz < 0) tz = (inner_z_ - local_point.z()) / local_dir.z(); + // Protect against outside the region. + if (tx < 0 || ty < 0 || tz < 0) + G4Exception("[BoxPointSampler]", "GetIntersect()", FatalException, + "Point outside region, projection fail!"); // The minimum of the tx, ty, tz gives the intersection point. G4double tmin = std::min({tx, ty, tz}); - return RotateAndTranslate(point + tmin * dir); + return RotateAndTranslate(local_point + tmin * local_dir); } } // end namespace nexus diff --git a/source/utils/BoxPointSampler.h b/source/utils/BoxPointSampler.h index 84e2349c52..ec5f5f21d1 100644 --- a/source/utils/BoxPointSampler.h +++ b/source/utils/BoxPointSampler.h @@ -39,6 +39,7 @@ namespace nexus { private: G4double GetLength(G4double origin, G4double max_length); G4ThreeVector RotateAndTranslate(G4ThreeVector position); + void InvertRotationAndTranslation(G4ThreeVector& vec, bool translate=true); private: G4double inner_x_, inner_y_, inner_z_; ///< Internal dimensions From 15a766ecb9d0fe25bb1706abb08abd725ebcdc7f Mon Sep 17 00:00:00 2001 From: Andrew Laing Date: Tue, 24 Nov 2020 18:32:27 +0100 Subject: [PATCH 07/25] Take into account arbitrary rotation/translation in CylinderSampler --- source/utils/CylinderPointSampler2020.cc | 54 ++++++++++++++++-------- source/utils/CylinderPointSampler2020.h | 1 + 2 files changed, 38 insertions(+), 17 deletions(-) diff --git a/source/utils/CylinderPointSampler2020.cc b/source/utils/CylinderPointSampler2020.cc index 01360031ce..aa1cd96d84 100644 --- a/source/utils/CylinderPointSampler2020.cc +++ b/source/utils/CylinderPointSampler2020.cc @@ -108,28 +108,40 @@ namespace nexus { CylinderPointSampler2020::GetIntersect(const G4ThreeVector& point, const G4ThreeVector& dir) { + // The point and direction should be in the lab frame + // Need to rotate into the Sampler frame to get the intersection + // Then rotate back into lab frame. + G4ThreeVector l_point = point; + InvertRotationAndTranslation(l_point); + G4ThreeVector l_dir = dir; + InvertRotationAndTranslation(l_dir, false); + G4cout << "Orig points: "<delta(), rotation_->axis()); + } } // end namespace nexus diff --git a/source/utils/CylinderPointSampler2020.h b/source/utils/CylinderPointSampler2020.h index d7e116aac5..cc7087379e 100644 --- a/source/utils/CylinderPointSampler2020.h +++ b/source/utils/CylinderPointSampler2020.h @@ -50,6 +50,7 @@ namespace nexus { G4double GetPhi(); G4double GetLength(G4double halfLength); G4ThreeVector RotateAndTranslate(G4ThreeVector position); + void InvertRotationAndTranslation(G4ThreeVector& vec, bool translate=true); private: G4double minRad_, maxRad_, halfLength_; // Solid Dimensions From 80a2dff1240c04208242536e1f5afc07483b1207 Mon Sep 17 00:00:00 2001 From: Andrew Laing Date: Tue, 24 Nov 2020 18:32:50 +0100 Subject: [PATCH 08/25] Cosmetics --- source/generators/MuonAngleGenerator.cc | 5 ++--- source/generators/MuonAngleGenerator.h | 1 - 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/source/generators/MuonAngleGenerator.cc b/source/generators/MuonAngleGenerator.cc index 6dacde506b..0aa4c4e36b 100644 --- a/source/generators/MuonAngleGenerator.cc +++ b/source/generators/MuonAngleGenerator.cc @@ -19,7 +19,6 @@ #include #include #include -#include "MuonsPointSampler.h" #include "AddUserInfoToPV.h" #include @@ -34,7 +33,7 @@ using namespace CLHEP; MuonAngleGenerator::MuonAngleGenerator(): G4VPrimaryGenerator(), msg_(0), particle_definition_(0), angular_generation_(true), rPhi_(NULL), energy_min_(0.), - energy_max_(0.), gen_rad_(223.33*cm), distribution_(0), geom_(0)//, geom_solid_(0) + energy_max_(0.), gen_rad_(223.33*cm), distribution_(0), geom_(0) { msg_ = new G4GenericMessenger(this, "/Generator/MuonAngleGenerator/", "Control commands of muongenerator."); @@ -188,7 +187,7 @@ void MuonAngleGenerator::GetDirection(G4ThreeVector& dir) G4ThreeVector MuonAngleGenerator::ProjectToVertex(const G4ThreeVector& dir) { - // Postion in disc (need to sort size, member function) + // Postion in disc G4double rad = gen_rad_ * std::sqrt(G4UniformRand()); G4double ang = 2 * G4UniformRand() * pi; diff --git a/source/generators/MuonAngleGenerator.h b/source/generators/MuonAngleGenerator.h index cb48f229a0..dd1394281e 100644 --- a/source/generators/MuonAngleGenerator.h +++ b/source/generators/MuonAngleGenerator.h @@ -16,7 +16,6 @@ class G4GenericMessenger; class G4Event; class G4ParticleDefinition; -//class G4VSolid; class TH2F; From 52bc173fddc34af6ca14055e5e71bd04163cc48b Mon Sep 17 00:00:00 2001 From: Andrew Laing Date: Thu, 26 Nov 2020 15:16:35 +0100 Subject: [PATCH 09/25] Add tests for intersection with samplers --- source/tests/utils/BoxPointSamplerTests.cc | 59 ++++++++++++++++++++++ source/utils/CylinderPointSampler2020.cc | 34 ++++++------- 2 files changed, 75 insertions(+), 18 deletions(-) diff --git a/source/tests/utils/BoxPointSamplerTests.cc b/source/tests/utils/BoxPointSamplerTests.cc index 8d6e6d7363..ec1f5c2412 100644 --- a/source/tests/utils/BoxPointSamplerTests.cc +++ b/source/tests/utils/BoxPointSamplerTests.cc @@ -50,3 +50,62 @@ TEST_CASE("BoxPointSampler") { } } + + +TEST_CASE("Expected intersect") { + // This test checks that the GetIntersect method + // of BoxPointSampler produces the expected + // intersect points for the simple cases of + // a central point with directions parallel to + // the faces of the box. + auto inner_dim = 100; + auto thickness = 10; + auto sampler = nexus::BoxPointSampler(inner_dim, inner_dim, inner_dim, thickness); + auto origin = G4ThreeVector(0., 0., 0.); + + G4double xdir[] = {1., 0., 0., -1., 0., 0.}; + G4double ydir[] = {0., 1., 0., 0., -1., 0.}; + G4double zdir[] = {0., 0., 1., 0., 0., -1.}; + for (G4int i=0; i<6; ++i){ + auto dir = G4ThreeVector(xdir[i], ydir[i], zdir[i]); + auto intersect = sampler.GetIntersect(origin, dir); + + REQUIRE(intersect.x() == Approx(inner_dim * dir.x())); + REQUIRE(intersect.y() == Approx(inner_dim * dir.y())); + REQUIRE(intersect.z() == Approx(inner_dim * dir.z())); + } + +} + +TEST_CASE("Box Arbitrary valid intersect") { + // This test checks that for a sampler of arbitrary + // position and arbitrary rotation a valid intersect + // point is found for a randomly generated ray. + auto inner_dim = 100; + auto thickness = 10; + + auto origin = G4ThreeVector(inner_dim * G4UniformRand(), + inner_dim * G4UniformRand(), + inner_dim * G4UniformRand()); + auto rotation = new G4RotationMatrix(); + rotation->rotateX(CLHEP::twopi * G4UniformRand()); + rotation->rotateY(CLHEP::twopi * G4UniformRand()); + rotation->rotateZ(CLHEP::twopi * G4UniformRand()); + + auto sampler = nexus::BoxPointSampler(inner_dim, inner_dim, inner_dim, thickness, + origin, rotation); + + auto point = G4ThreeVector(inner_dim * G4UniformRand(), + inner_dim * G4UniformRand(), + inner_dim * G4UniformRand()); + auto dir = G4ThreeVector(G4UniformRand(), + G4UniformRand(), + G4UniformRand()).unit(); + + auto intersect = sampler.GetIntersect(point, dir); + + // Check that the new ray passes through the original point. + auto origin_point = point - intersect; + REQUIRE(std::abs(origin_point.dot(dir)) == Approx(origin_point.mag())); + +} diff --git a/source/utils/CylinderPointSampler2020.cc b/source/utils/CylinderPointSampler2020.cc index aa1cd96d84..88af6d8ebd 100644 --- a/source/utils/CylinderPointSampler2020.cc +++ b/source/utils/CylinderPointSampler2020.cc @@ -115,8 +115,6 @@ namespace nexus { InvertRotationAndTranslation(l_point); G4ThreeVector l_dir = dir; InvertRotationAndTranslation(l_dir, false); - G4cout << "Orig points: "< Date: Wed, 9 Dec 2020 14:53:22 +0100 Subject: [PATCH 10/25] Add tests for CylinderPointSampler2020 --- .../tests/utils/CylinderPointSamplerTests.cc | 76 +++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100644 source/tests/utils/CylinderPointSamplerTests.cc diff --git a/source/tests/utils/CylinderPointSamplerTests.cc b/source/tests/utils/CylinderPointSamplerTests.cc new file mode 100644 index 0000000000..c58ddecaaf --- /dev/null +++ b/source/tests/utils/CylinderPointSamplerTests.cc @@ -0,0 +1,76 @@ +#include +#include + +#include + +TEST_CASE("Barrel and Caps intersection") { + // Tests that rays from the origin parallel + // to the axes are found to intersect with the + // Cylinder sampler at the expected points. + auto minRad = 200; + auto maxRad = 220; + auto halfLen = 500; + auto sampler = nexus::CylinderPointSampler2020(minRad, maxRad, halfLen); + auto origin = G4ThreeVector(0., 0., 0.); + + // Check caps. + auto intersect = sampler.GetIntersect(origin, G4ThreeVector(0., 0., 1.)); + REQUIRE(intersect.x() == Approx(0.)); + REQUIRE(intersect.y() == Approx(0.)); + REQUIRE(intersect.z() == Approx(halfLen)); + intersect = sampler.GetIntersect(origin, G4ThreeVector(0., 0., -1.)); + REQUIRE(intersect.x() == Approx(0.)); + REQUIRE(intersect.y() == Approx(0.)); + REQUIRE(intersect.z() == Approx(-halfLen)); + + // Now the barrel. + intersect = sampler.GetIntersect(origin, G4ThreeVector(1., 0., 0.)); + REQUIRE(intersect.x() == Approx(minRad)); + REQUIRE(intersect.y() == Approx(0.)); + REQUIRE(intersect.z() == Approx(0.)); + intersect = sampler.GetIntersect(origin, G4ThreeVector(-1., 0., 0.)); + REQUIRE(intersect.x() == Approx(-minRad)); + REQUIRE(intersect.y() == Approx(0.)); + REQUIRE(intersect.z() == Approx(0.)); + intersect = sampler.GetIntersect(origin, G4ThreeVector(0., 1., 0.)); + REQUIRE(intersect.x() == Approx(0.)); + REQUIRE(intersect.y() == Approx(minRad)); + REQUIRE(intersect.z() == Approx(0.)); + intersect = sampler.GetIntersect(origin, G4ThreeVector(0., -1., 0.)); + REQUIRE(intersect.x() == Approx(0.)); + REQUIRE(intersect.y() == Approx(-minRad)); + REQUIRE(intersect.z() == Approx(0.)); + +} + +TEST_CASE("Cylinder Arbitrary valid intersect") { + // This test checks that for a sampler of arbitrary + // position and arbitrary rotation a valid intersect + // point is found for a randomly generated ray. + auto minRad = 200; + auto maxRad = 220; + auto halfLen = 500; + auto origin = G4ThreeVector(maxRad * G4UniformRand(), + maxRad * G4UniformRand(), + halfLen * G4UniformRand()); + auto rotation = new G4RotationMatrix(); + rotation->rotateX(CLHEP::twopi * G4UniformRand()); + rotation->rotateY(CLHEP::twopi * G4UniformRand()); + rotation->rotateZ(CLHEP::twopi * G4UniformRand()); + auto sampler = nexus::CylinderPointSampler2020(minRad, maxRad, halfLen, + 0., CLHEP::twopi, + rotation, origin); + auto point = G4ThreeVector(minRad * G4UniformRand(), + minRad * G4UniformRand(), + minRad * G4UniformRand()); + auto dir = G4ThreeVector(G4UniformRand(), + G4UniformRand(), + G4UniformRand()).unit(); + + auto intersect = sampler.GetIntersect(point, dir); + + // Check that the new ray passes through the original point. + auto origin_point = point - intersect; + REQUIRE(std::abs(origin_point.dot(dir)) == Approx(origin_point.mag())); + +} From a2db2e46d8a7a0d47c7bb3765fb8d1d43dcda1d0 Mon Sep 17 00:00:00 2001 From: Andrew Laing Date: Wed, 13 Jan 2021 16:47:30 +0100 Subject: [PATCH 11/25] Add ProjectToRegion funciton to LSCHallA class --- source/geometries/LSCHallA.cc | 19 +++++++++++++++++++ source/geometries/LSCHallA.h | 6 ++++++ 2 files changed, 25 insertions(+) diff --git a/source/geometries/LSCHallA.cc b/source/geometries/LSCHallA.cc index c9d3d5df35..3f1b9a35be 100644 --- a/source/geometries/LSCHallA.cc +++ b/source/geometries/LSCHallA.cc @@ -119,4 +119,23 @@ namespace nexus { return vertex; } + G4ThreeVector LSCHallA::ProjectToRegion(const G4String& region, + const G4ThreeVector& point, + const G4ThreeVector& dir) const + { + // Project backwards along dir from point to find the first intersection + // with region. + G4ThreeVector vertex(0., 0., 0.); + if (region == "HALLA_INNER") + return hallA_vertex_gen_->GetIntersect(point, dir); + else if (region == "HALLA_OUTER") + return hallA_outer_gen_->GetIntersect(point, dir); + else { + G4Exception("[LSCHallA]", "ProjectToRegion()", FatalException, + "Unknown vertex generation region!"); + } + + return vertex; + } + } diff --git a/source/geometries/LSCHallA.h b/source/geometries/LSCHallA.h index fc1a10c12d..0017f44809 100644 --- a/source/geometries/LSCHallA.h +++ b/source/geometries/LSCHallA.h @@ -30,6 +30,12 @@ namespace nexus { /// Generate a vertex within a given region of the geometry G4ThreeVector GenerateVertex(const G4String& region) const; + /// Returns a point within a region projecting from a + /// given point backwards along a line. + G4ThreeVector ProjectToRegion(const G4String& region, + const G4ThreeVector& point, + const G4ThreeVector& dir) const; + /// Builder void Construct(); From 2ed112e4d4b5474cea3bcdb8690dbb6b09a79f7e Mon Sep 17 00:00:00 2001 From: Andrew Laing Date: Wed, 13 Jan 2021 16:48:53 +0100 Subject: [PATCH 12/25] Allow HallA regions in ProjectToRegion for NextNew and Next100 --- source/geometries/Next100.cc | 6 ++++++ source/geometries/NextNew.cc | 6 ++++++ 2 files changed, 12 insertions(+) diff --git a/source/geometries/Next100.cc b/source/geometries/Next100.cc index 446c5ed1dd..348645998d 100644 --- a/source/geometries/Next100.cc +++ b/source/geometries/Next100.cc @@ -282,6 +282,12 @@ namespace nexus { if (region == "EXTERNAL"){ return shielding_->ProjectToRegion(region, point, dir); } + else if ((region == "HALLA_OUTER") || (region == "HALLA_INNER")){ + if (!lab_walls_) + G4Exception("[Next100]", "ProjectToRegion()", FatalException, + "To project to this region you need lab_walls == true!"); + return hallA_walls_->ProjectToRegion(region, point, dir); + } else { G4Exception("[Next100]", "ProjectToRegion()", FatalException, "Unknown vertex generation region!"); diff --git a/source/geometries/NextNew.cc b/source/geometries/NextNew.cc index e33e3df730..b922bb3f36 100644 --- a/source/geometries/NextNew.cc +++ b/source/geometries/NextNew.cc @@ -609,6 +609,12 @@ namespace nexus { if (region == "EXTERNAL"){ return shielding_->ProjectToRegion(region, point, dir); } + else if ((region == "HALLA_OUTER") || (region == "HALLA_INNER")){ + if (!lab_walls_) + G4Exception("[NextNew]", "ProjectToRegion()", FatalException, + "To project to this region you need lab_walls == true!"); + return hallA_walls_->ProjectToRegion(region, point, dir); + } else { G4Exception("[NextNew]", "ProjectToRegion()", FatalException, "Unknown vertex generation region!"); From 47abb2e69827d2fcf8f22e45f1158ed409d37d36 Mon Sep 17 00:00:00 2001 From: Andrew Laing Date: Wed, 13 Jan 2021 16:52:31 +0100 Subject: [PATCH 13/25] Fix explicative comments in ProjectToRegion functions --- source/geometries/LSCHallA.cc | 2 +- source/geometries/Next100.cc | 2 +- source/geometries/Next100Shielding.cc | 2 +- source/geometries/NextNew.cc | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/source/geometries/LSCHallA.cc b/source/geometries/LSCHallA.cc index 3f1b9a35be..fd7672b220 100644 --- a/source/geometries/LSCHallA.cc +++ b/source/geometries/LSCHallA.cc @@ -123,7 +123,7 @@ namespace nexus { const G4ThreeVector& point, const G4ThreeVector& dir) const { - // Project backwards along dir from point to find the first intersection + // Project along dir from point to find the first intersection // with region. G4ThreeVector vertex(0., 0., 0.); if (region == "HALLA_INNER") diff --git a/source/geometries/Next100.cc b/source/geometries/Next100.cc index 348645998d..8af475070c 100644 --- a/source/geometries/Next100.cc +++ b/source/geometries/Next100.cc @@ -276,7 +276,7 @@ namespace nexus { const G4ThreeVector& point, const G4ThreeVector& dir) const { - // Project backwards along dir from point to find the first intersection + // Project along dir from point to find the first intersection // with region. G4ThreeVector vertex(0., 0., 0.); if (region == "EXTERNAL"){ diff --git a/source/geometries/Next100Shielding.cc b/source/geometries/Next100Shielding.cc index 18c87de068..29deb5d168 100644 --- a/source/geometries/Next100Shielding.cc +++ b/source/geometries/Next100Shielding.cc @@ -800,7 +800,7 @@ namespace nexus { const G4ThreeVector& point, const G4ThreeVector& dir) const { - // Project backwards along dir from point to find the first intersection + // Project along dir from point to find the first intersection // with region. G4ThreeVector vertex(0., 0., 0.); if (region == "EXTERNAL"){ diff --git a/source/geometries/NextNew.cc b/source/geometries/NextNew.cc index b922bb3f36..7a9321eada 100644 --- a/source/geometries/NextNew.cc +++ b/source/geometries/NextNew.cc @@ -603,7 +603,7 @@ namespace nexus { const G4ThreeVector& point, const G4ThreeVector& dir) const { - // Project backwards along dir from point to find the first intersection + // Project along dir from point to find the first intersection // with region. G4ThreeVector vertex(0., 0., 0.); if (region == "EXTERNAL"){ From 5dfe78637201b6dbda88ae7d6dcfe832f4242fcb Mon Sep 17 00:00:00 2001 From: Paola Ferrario Date: Thu, 28 Sep 2023 09:47:15 +0200 Subject: [PATCH 14/25] Clarify function for macro NEXT100_muons_lsc --- macros/NEXT100_muons.init.mac | 4 ++- macros/NEXT100_muons_lsc.config.mac | 44 +++++++++++++++++++++++++++++ 2 files changed, 47 insertions(+), 1 deletion(-) create mode 100644 macros/NEXT100_muons_lsc.config.mac diff --git a/macros/NEXT100_muons.init.mac b/macros/NEXT100_muons.init.mac index 6bda9da6a6..7e6cd645d2 100644 --- a/macros/NEXT100_muons.init.mac +++ b/macros/NEXT100_muons.init.mac @@ -1,7 +1,9 @@ ## ---------------------------------------------------------------------------- ## nexus | NEXT100_muons_lsc.init.mac ## -## Initialization macro to simulate muons in the NEXT-100 geometry. +## Initialization macro to simulate muons according to the distribution +## measured at LSC, in the NEXT-100 geometry launching from just outside +## the lead shielding. ## ## The NEXT Collaboration ## ---------------------------------------------------------------------------- diff --git a/macros/NEXT100_muons_lsc.config.mac b/macros/NEXT100_muons_lsc.config.mac new file mode 100644 index 0000000000..d109cc8cdb --- /dev/null +++ b/macros/NEXT100_muons_lsc.config.mac @@ -0,0 +1,44 @@ +## ---------------------------------------------------------------------------- +## nexus | NEXT100_muons_lsc.config.mac +## +## Configuration macro to simulate muons according to the distribution +## measured at LSC, in the NEXT-100 geometry launching from just outside +## the lead shielding. +## +## The NEXT Collaboration +## ---------------------------------------------------------------------------- + +### VERBOSITY +/control/verbose 0 +/run/verbose 0 +/event/verbose 0 +/tracking/verbose 0 + +### JOB CONTROL +/nexus/random_seed 17392 + +### GEOMETRY +/Geometry/Next100/pressure 15 bar +/Geometry/Next100/gas enrichedXe +/Geometry/Next100/elfield false + +### GENERATOR +/Generator/MuonAngleGenerator/region EXTERNAL +/Generator/MuonAngleGenerator/min_energy 200 GeV +/Generator/MuonAngleGenerator/max_energy 250 GeV +/Generator/MuonAngleGenerator/azimuth_rotation 150 deg +/Generator/MuonAngleGenerator/angle_file histos/MuonAnaAllRuns.root +/Generator/MuonAngleGenerator/angle_dist za + +### ACTIONS +/Actions/DefaultEventAction/energy_threshold 0.01 MeV + + +### PHYSICS (for fast simulation) +/PhysicsList/Nexus/clustering false +/PhysicsList/Nexus/drift false +/PhysicsList/Nexus/electroluminescence false + +### PERSISTENCY +/nexus/persistency/start_id 0 +/nexus/persistency/outputFile Next100Muons_lsc_example.next \ No newline at end of file From a4e1192a5ec9ddf303c324f642d4737d9b750654 Mon Sep 17 00:00:00 2001 From: Andrew Laing Date: Tue, 26 Jan 2021 18:28:57 +0100 Subject: [PATCH 15/25] Guarantee generated point inside volume in GetIntersect tests --- source/tests/utils/BoxPointSamplerTests.cc | 8 +++----- source/tests/utils/CylinderPointSamplerTests.cc | 9 ++++++--- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/source/tests/utils/BoxPointSamplerTests.cc b/source/tests/utils/BoxPointSamplerTests.cc index ec1f5c2412..7034f42113 100644 --- a/source/tests/utils/BoxPointSamplerTests.cc +++ b/source/tests/utils/BoxPointSamplerTests.cc @@ -92,12 +92,10 @@ TEST_CASE("Box Arbitrary valid intersect") { rotation->rotateY(CLHEP::twopi * G4UniformRand()); rotation->rotateZ(CLHEP::twopi * G4UniformRand()); - auto sampler = nexus::BoxPointSampler(inner_dim, inner_dim, inner_dim, thickness, - origin, rotation); + auto sampler = nexus::BoxPointSampler(inner_dim, inner_dim, inner_dim, + thickness, origin, rotation); - auto point = G4ThreeVector(inner_dim * G4UniformRand(), - inner_dim * G4UniformRand(), - inner_dim * G4UniformRand()); + auto point = sampler.GenerateVertex("INSIDE"); auto dir = G4ThreeVector(G4UniformRand(), G4UniformRand(), G4UniformRand()).unit(); diff --git a/source/tests/utils/CylinderPointSamplerTests.cc b/source/tests/utils/CylinderPointSamplerTests.cc index c58ddecaaf..d9d480340f 100644 --- a/source/tests/utils/CylinderPointSamplerTests.cc +++ b/source/tests/utils/CylinderPointSamplerTests.cc @@ -60,9 +60,12 @@ TEST_CASE("Cylinder Arbitrary valid intersect") { auto sampler = nexus::CylinderPointSampler2020(minRad, maxRad, halfLen, 0., CLHEP::twopi, rotation, origin); - auto point = G4ThreeVector(minRad * G4UniformRand(), - minRad * G4UniformRand(), - minRad * G4UniformRand()); + auto point = G4ThreeVector(minRad * (2 * G4UniformRand() - 1), + minRad * (2 * G4UniformRand() - 1), + halfLen * (2 * G4UniformRand() - 1)); + // Rotate and translate so inside Sampler Volume. + point *= *rotation; + point += origin; auto dir = G4ThreeVector(G4UniformRand(), G4UniformRand(), G4UniformRand()).unit(); From afdc4640bbeacde0e64a7a09e55bad0098b64850 Mon Sep 17 00:00:00 2001 From: Andrew Laing Date: Tue, 26 Jan 2021 18:42:42 +0100 Subject: [PATCH 16/25] Add additional explanation of generation method --- source/generators/MuonAngleGenerator.cc | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/source/generators/MuonAngleGenerator.cc b/source/generators/MuonAngleGenerator.cc index 0aa4c4e36b..4a554fa007 100644 --- a/source/generators/MuonAngleGenerator.cc +++ b/source/generators/MuonAngleGenerator.cc @@ -187,12 +187,26 @@ void MuonAngleGenerator::GetDirection(G4ThreeVector& dir) G4ThreeVector MuonAngleGenerator::ProjectToVertex(const G4ThreeVector& dir) { + ///////////////////////////////////////////////////////////////////////// + // This method of vertex generation decides the starting vertex + // of the primary particle in three steps: + // 1) A random point is generated on a disc centred on the main + // volumes. + // 2) This disc is rotated so that it is perpendicular to the + // direction vector which is the only argument of the function. + // The new point is a point through which the vector passes. + // 3) The vertex is found by projecting backwards along the direction + // vector from that point to find the point where the ray + // (point - t*dir) intersects with the region configured as the + // starting point for all vertices. + ///////////////////////////////////////////////////////////////////////// // Postion in disc G4double rad = gen_rad_ * std::sqrt(G4UniformRand()); G4double ang = 2 * G4UniformRand() * pi; // Continue assuming that Y is vertical and z drift, - // valid for NEW and NEXT-100 (at least) + // valid for NEW and NEXT-100 (at least). + // Rotate the disc (origin-point vector) to be perpendicular to dir. G4ThreeVector point(rad * std::cos(ang), 0., rad * std::sin(ang)); point.rotate(pi / 2 - dir.angle(point), dir.cross(point)); From d1ae9417e43ef9baed6d6d98d091354fb9575982 Mon Sep 17 00:00:00 2001 From: Paola Ferrario Date: Thu, 28 Sep 2023 09:48:59 +0200 Subject: [PATCH 17/25] Eliminate obsolete macro --- macros/NEXT100_muons_lsc.config.mac | 44 ----------------------------- 1 file changed, 44 deletions(-) delete mode 100644 macros/NEXT100_muons_lsc.config.mac diff --git a/macros/NEXT100_muons_lsc.config.mac b/macros/NEXT100_muons_lsc.config.mac deleted file mode 100644 index d109cc8cdb..0000000000 --- a/macros/NEXT100_muons_lsc.config.mac +++ /dev/null @@ -1,44 +0,0 @@ -## ---------------------------------------------------------------------------- -## nexus | NEXT100_muons_lsc.config.mac -## -## Configuration macro to simulate muons according to the distribution -## measured at LSC, in the NEXT-100 geometry launching from just outside -## the lead shielding. -## -## The NEXT Collaboration -## ---------------------------------------------------------------------------- - -### VERBOSITY -/control/verbose 0 -/run/verbose 0 -/event/verbose 0 -/tracking/verbose 0 - -### JOB CONTROL -/nexus/random_seed 17392 - -### GEOMETRY -/Geometry/Next100/pressure 15 bar -/Geometry/Next100/gas enrichedXe -/Geometry/Next100/elfield false - -### GENERATOR -/Generator/MuonAngleGenerator/region EXTERNAL -/Generator/MuonAngleGenerator/min_energy 200 GeV -/Generator/MuonAngleGenerator/max_energy 250 GeV -/Generator/MuonAngleGenerator/azimuth_rotation 150 deg -/Generator/MuonAngleGenerator/angle_file histos/MuonAnaAllRuns.root -/Generator/MuonAngleGenerator/angle_dist za - -### ACTIONS -/Actions/DefaultEventAction/energy_threshold 0.01 MeV - - -### PHYSICS (for fast simulation) -/PhysicsList/Nexus/clustering false -/PhysicsList/Nexus/drift false -/PhysicsList/Nexus/electroluminescence false - -### PERSISTENCY -/nexus/persistency/start_id 0 -/nexus/persistency/outputFile Next100Muons_lsc_example.next \ No newline at end of file From e9c063123905b42b65f8c06003d2b4fe97201b9b Mon Sep 17 00:00:00 2001 From: Paola Ferrario Date: Thu, 28 Sep 2023 09:57:49 +0200 Subject: [PATCH 18/25] Adapt new functions to new base class --- source/geometries/GeometryBase.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/geometries/GeometryBase.h b/source/geometries/GeometryBase.h index 6a70b4c952..cf07b34524 100644 --- a/source/geometries/GeometryBase.h +++ b/source/geometries/GeometryBase.h @@ -108,12 +108,12 @@ namespace nexus { inline G4ThreeVector GeometryBase::GenerateVertex(const G4String&) const { return G4ThreeVector(0., 0., 0.); } - inline G4ThreeVector BaseGeometry::ProjectToRegion(const G4String&, + inline G4ThreeVector GeometryBase::ProjectToRegion(const G4String&, const G4ThreeVector&, const G4ThreeVector&) const { return G4ThreeVector(0., 0., 0.); } - inline void BaseGeometry::SetSpan(G4double s) { span_ = s; } + inline void GeometryBase::SetSpan(G4double s) { span_ = s; } inline G4double GeometryBase::GetSpan() { return span_; } From 13ea4a02bb47eb5b34508861447aa30a95c54b93 Mon Sep 17 00:00:00 2001 From: Paola Ferrario Date: Thu, 28 Sep 2023 09:58:47 +0200 Subject: [PATCH 19/25] Remove obsolete class --- source/generators/MuonAngleGenerator.cc | 215 ------------------------ source/generators/MuonAngleGenerator.h | 81 --------- 2 files changed, 296 deletions(-) delete mode 100644 source/generators/MuonAngleGenerator.cc delete mode 100644 source/generators/MuonAngleGenerator.h diff --git a/source/generators/MuonAngleGenerator.cc b/source/generators/MuonAngleGenerator.cc deleted file mode 100644 index 4a554fa007..0000000000 --- a/source/generators/MuonAngleGenerator.cc +++ /dev/null @@ -1,215 +0,0 @@ -// ---------------------------------------------------------------------------- -// nexus | MuonAngleGenerator.cc -// -// This class is the primary generator of muons following an angular -// distribution measured in the LSC. -// -// The NEXT Collaboration -// ---------------------------------------------------------------------------- - -#include "MuonAngleGenerator.h" -#include "DetectorConstruction.h" -#include "BaseGeometry.h" -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include "AddUserInfoToPV.h" - -#include -#include "TFile.h" -#include "TH2F.h" -#include "CLHEP/Units/SystemOfUnits.h" - -using namespace nexus; -using namespace CLHEP; - - -MuonAngleGenerator::MuonAngleGenerator(): - G4VPrimaryGenerator(), msg_(0), particle_definition_(0), - angular_generation_(true), rPhi_(NULL), energy_min_(0.), - energy_max_(0.), gen_rad_(223.33*cm), distribution_(0), geom_(0) -{ - msg_ = new G4GenericMessenger(this, "/Generator/MuonAngleGenerator/", - "Control commands of muongenerator."); - - G4GenericMessenger::Command& min_energy = - msg_->DeclareProperty("min_energy", energy_min_, "Set minimum kinetic energy of the particle."); - min_energy.SetUnitCategory("Energy"); - min_energy.SetParameterName("min_energy", false); - min_energy.SetRange("min_energy>0."); - - G4GenericMessenger::Command& max_energy = - msg_->DeclareProperty("max_energy", energy_max_, "Set maximum kinetic energy of the particle"); - max_energy.SetUnitCategory("Energy"); - max_energy.SetParameterName("max_energy", false); - max_energy.SetRange("max_energy>0."); - - // defaults to length of the NEXT100 shielding diagonal if not present. - G4GenericMessenger::Command& generation_radius = - msg_->DeclareProperty("gen_rad", gen_rad_, "Set radius for generation disc"); - generation_radius.SetUnitCategory("Length"); - generation_radius.SetParameterName("gen_rad", false); - generation_radius.SetRange("gen_rad>0."); - - msg_->DeclareProperty("region", region_, - "Set the region of the geometry where the vertex will be generated."); - - msg_->DeclareProperty("angles_on", angular_generation_, - "Distribute muon directions according to file?"); - - msg_->DeclareProperty("angle_file", ang_file_, - "Name of the file containing angular distribution."); - msg_->DeclareProperty("angle_dist", dist_name_, - "Name of the angular distribution histogram."); - - G4GenericMessenger::Command& rotation = - msg_->DeclareProperty("azimuth_rotation", axis_rotation_, - "Angle between north and nexus z in anticlockwise"); - rotation.SetUnitCategory("Angle"); - rotation.SetParameterName("azimuth", false); - rotation.SetRange("azimuth>0."); - - DetectorConstruction* detconst = (DetectorConstruction*) G4RunManager::GetRunManager()->GetUserDetectorConstruction(); - geom_ = detconst->GetGeometry(); - -} - - -MuonAngleGenerator::~MuonAngleGenerator() -{ - - delete msg_; -} - - -void MuonAngleGenerator::SetupAngles() -{ - // Rotation from the axes used in file. - // Rotates anticlockwise about Y. - rPhi_ = new G4RotationMatrix(); - rPhi_->rotateY(-axis_rotation_); - - // Get the Angular distribution from file. - TFile angle_file(ang_file_); - angle_file.GetObject(dist_name_, distribution_); - distribution_->SetDirectory(0); - angle_file.Close(); - -} - - -void MuonAngleGenerator::GeneratePrimaryVertex(G4Event* event) -{ - - if (angular_generation_ && rPhi_ == NULL) - SetupAngles(); - - particle_definition_ = - G4ParticleTable::GetParticleTable()->FindParticle(MuonCharge()); - if (!particle_definition_) - G4Exception("[MuonAngleGenerator]", "SetParticleDefinition()", - FatalException, " can not create a muon "); - - // Generate uniform random energy in [E_min, E_max] - G4double kinetic_energy = RandomEnergy(); - - // Particle propierties - G4double mass = particle_definition_->GetPDGMass(); - G4double energy = kinetic_energy + mass; - G4double pmod = std::sqrt(energy*energy - mass*mass); - - G4ThreeVector p_dir(0., -1., 0.); - if (angular_generation_) - GetDirection(p_dir); - G4ThreeVector position = ProjectToVertex(p_dir); - G4double px = pmod * p_dir.x(); - G4double py = pmod * p_dir.y(); - G4double pz = pmod * p_dir.z(); - - // Particle generated at start-of-event - G4double time = 0.; - // Create a new vertex - G4PrimaryVertex* vertex = new G4PrimaryVertex(position, time); - - // Create the new primary particle and set it some properties - G4PrimaryParticle* particle = - new G4PrimaryParticle(particle_definition_, px, py, pz); - - // Add particle to the vertex and this to the event - vertex->SetPrimary(particle); - event->AddPrimaryVertex(vertex); -} - - -G4double MuonAngleGenerator::RandomEnergy() const -{ - if (energy_max_ == energy_min_) - return energy_min_; - else - return (G4UniformRand()*(energy_max_ - energy_min_) + energy_min_); -} - - -G4String MuonAngleGenerator::MuonCharge() const -{ - G4double rndCh = 2.3 *G4UniformRand(); //From PDG cosmic muons mu+/mu- = 1.3 - if (rndCh <1.3) - return "mu+"; - else - return "mu-"; -} - - -void MuonAngleGenerator::GetDirection(G4ThreeVector& dir) -{ - // GetAngles from file?? Azimuth defined anticlockwise - // From north - G4double zenith = 0.; - G4double azimuth = 0.; - distribution_->GetRandom2(azimuth, zenith); - // !! Current distribution in units of pi - zenith *= pi; - azimuth *= pi; - - dir.setX(sin(zenith) * sin(azimuth)); - dir.setY(-cos(zenith)); - dir.setZ(-sin(zenith) * cos(azimuth)); - - dir *= *rPhi_; -} - - -G4ThreeVector MuonAngleGenerator::ProjectToVertex(const G4ThreeVector& dir) -{ - ///////////////////////////////////////////////////////////////////////// - // This method of vertex generation decides the starting vertex - // of the primary particle in three steps: - // 1) A random point is generated on a disc centred on the main - // volumes. - // 2) This disc is rotated so that it is perpendicular to the - // direction vector which is the only argument of the function. - // The new point is a point through which the vector passes. - // 3) The vertex is found by projecting backwards along the direction - // vector from that point to find the point where the ray - // (point - t*dir) intersects with the region configured as the - // starting point for all vertices. - ///////////////////////////////////////////////////////////////////////// - // Postion in disc - G4double rad = gen_rad_ * std::sqrt(G4UniformRand()); - G4double ang = 2 * G4UniformRand() * pi; - - // Continue assuming that Y is vertical and z drift, - // valid for NEW and NEXT-100 (at least). - // Rotate the disc (origin-point vector) to be perpendicular to dir. - G4ThreeVector point(rad * std::cos(ang), 0., rad * std::sin(ang)); - point.rotate(pi / 2 - dir.angle(point), dir.cross(point)); - - // Now project back to the requested region intersection. - return geom_->ProjectToRegion(region_, point, -dir); -} diff --git a/source/generators/MuonAngleGenerator.h b/source/generators/MuonAngleGenerator.h deleted file mode 100644 index dd1394281e..0000000000 --- a/source/generators/MuonAngleGenerator.h +++ /dev/null @@ -1,81 +0,0 @@ -// ---------------------------------------------------------------------------- -// nexus | MuonAngleGenerator.h -// -// This class is the primary generator of muons following an angular -// distribution measured in the LSC. -// -// The NEXT Collaboration -// ---------------------------------------------------------------------------- - -#ifndef MUON_ANGLE_GENERATOR_H -#define MUON_ANGLE_GENERATOR_H - -#include -#include - -class G4GenericMessenger; -class G4Event; -class G4ParticleDefinition; - -class TH2F; - - -namespace nexus { - - class BaseGeometry; - - class MuonAngleGenerator: public G4VPrimaryGenerator - { - public: - /// Constructor - MuonAngleGenerator(); - /// Destructor - ~MuonAngleGenerator(); - - /// This method is invoked at the beginning of the event. It sets - /// a primary vertex (that is, a particle in a given position and time) - /// in the event. - void GeneratePrimaryVertex(G4Event*); - - private: - - // Sets the rotation angle and the spectra to - // be read for angle generation as well as - // setting the overlap volume for filtering. - void SetupAngles(); - - /// Generate a random kinetic energy with flat probability in - // the interval [energy_min, energy_max]. - G4double RandomEnergy() const; - G4String MuonCharge() const; - - void GetDirection(G4ThreeVector& dir); - - G4ThreeVector ProjectToVertex(const G4ThreeVector& dir); - - private: - G4GenericMessenger* msg_; - - G4ParticleDefinition* particle_definition_; - - G4bool angular_generation_; ///< Distribution or all downwards - G4double axis_rotation_; ///< Angle between North and +z - G4RotationMatrix *rPhi_; ///< Rotation to adjust axes - - G4double energy_min_; ///< Minimum kinetic energy - G4double energy_max_; ///< Maximum kinetic energy - - G4String region_; ///< Name of generator region - G4String ang_file_; ///< Name of file with distributions - G4String dist_name_; ///< Name of distribution in file - G4double gen_rad_; ///< Radius of disc for generation - - TH2F * distribution_; ///< Anglular distribution - - const BaseGeometry* geom_; ///< Pointer to the detector geometry - - }; - -} // end namespace nexus - -#endif From 3d975ee162b2a9c467f6848366015d4634c96d22 Mon Sep 17 00:00:00 2001 From: Paola Ferrario Date: Thu, 28 Sep 2023 10:07:22 +0200 Subject: [PATCH 20/25] Use offset from maximum diagonal --- source/geometries/Next100Shielding.cc | 4 ---- 1 file changed, 4 deletions(-) diff --git a/source/geometries/Next100Shielding.cc b/source/geometries/Next100Shielding.cc index 29deb5d168..189b04a210 100644 --- a/source/geometries/Next100Shielding.cc +++ b/source/geometries/Next100Shielding.cc @@ -384,10 +384,6 @@ namespace nexus { inner_air_gen_ = new BoxPointSampler(shield_x_, shield_y_, shield_z_, 0, G4ThreeVector(0., 0., 0.), 0); - G4double ext_offset = 1. * cm; - external_gen_ = - new BoxPointSampler(lead_x_ + ext_offset, lead_y_ + ext_offset, lead_z_ + ext_offset, 1. * mm, - G4ThreeVector(0., 0., 0.), 0); // STEEL STRUCTURE GENERATORS lat_roof_gen_ = From 6e6234d069d465641343603f8122e2bf2629d530 Mon Sep 17 00:00:00 2001 From: Paola Ferrario Date: Thu, 28 Sep 2023 10:14:40 +0200 Subject: [PATCH 21/25] Merge optimized code to generate muons --- source/generators/MuonGenerator.cc | 50 +++++++++++++++++++++++++++--- source/generators/MuonGenerator.h | 4 +++ 2 files changed, 49 insertions(+), 5 deletions(-) diff --git a/source/generators/MuonGenerator.cc b/source/generators/MuonGenerator.cc index 7dc5c156a0..158f474e1d 100644 --- a/source/generators/MuonGenerator.cc +++ b/source/generators/MuonGenerator.cc @@ -34,7 +34,8 @@ REGISTER_CLASS(MuonGenerator, G4VPrimaryGenerator) MuonGenerator::MuonGenerator(): G4VPrimaryGenerator(), msg_(0), particle_definition_(0), use_lsc_dist_(true), axis_rotation_(150), rPhi_(NULL), user_dir_{}, energy_min_(0.), - energy_max_(0.), dist_name_("za"), bInitialize_(false), geom_(0), geom_solid_(0) + energy_max_(0.), dist_name_("za"), bInitialize_(false), geom_(0), geom_solid_(0), + gen_rad_(223.33*cm) { msg_ = new G4GenericMessenger(this, "/Generator/MuonGenerator/", "Control commands of muongenerator."); @@ -71,7 +72,15 @@ MuonGenerator::MuonGenerator(): rotation.SetParameterName("azimuth", false); rotation.SetRange("azimuth>0."); - DetectorConstruction* detconst = (DetectorConstruction*) G4RunManager::GetRunManager()->GetUserDetectorConstruction(); + // defaults to length of the NEXT100 shielding diagonal if not present. + G4GenericMessenger::Command& generation_radius = + msg_->DeclareProperty("gen_rad", gen_rad_, "Set radius for generation disc"); + generation_radius.SetUnitCategory("Length"); + generation_radius.SetParameterName("gen_rad", false); + generation_radius.SetRange("gen_rad>0."); + + DetectorConstruction* detconst = + (DetectorConstruction*) G4RunManager::GetRunManager()->GetUserDetectorConstruction(); geom_ = detconst->GetGeometry(); } @@ -201,7 +210,7 @@ void MuonGenerator::GeneratePrimaryVertex(G4Event* event) // Particle properties G4double mass = particle_definition_->GetPDGMass(); G4double energy = kinetic_energy + mass; - G4ThreeVector position = geom_->GenerateVertex(region_); + //G4ThreeVector position = geom_->GenerateVertex(region_); // Set default momentum and angular variables G4ThreeVector p_dir; @@ -211,7 +220,7 @@ void MuonGenerator::GeneratePrimaryVertex(G4Event* event) // Momentum, zenith, azimuth (and energy) from angular distribution file if (use_lsc_dist_){ GetDirection(p_dir, zenith, azimuth, energy, kinetic_energy, mass); - position = geom_->GenerateVertex(region_); + // position = geom_->GenerateVertex(region_); } else { @@ -236,9 +245,10 @@ void MuonGenerator::GeneratePrimaryVertex(G4Event* event) p_dir *= *rPhi_; } - } + G4ThreeVector position = ProjectToVertex(p_dir); + G4double pmod = std::sqrt(energy*energy - mass*mass); G4double px = pmod * p_dir.x(); G4double py = pmod * p_dir.y(); @@ -327,6 +337,36 @@ void MuonGenerator::GetDirection(G4ThreeVector& dir, G4double& zenith, G4double& } +G4ThreeVector MuonGenerator::ProjectToVertex(const G4ThreeVector& dir) +{ + ///////////////////////////////////////////////////////////////////////// + // This method of vertex generation decides the starting vertex + // of the primary particle in three steps: + // 1) A random point is generated on a disc centred on the main + // volumes. + // 2) This disc is rotated so that it is perpendicular to the + // direction vector which is the only argument of the function. + // The new point is a point through which the vector passes. + // 3) The vertex is found by projecting backwards along the direction + // vector from that point to find the point where the ray + // (point - t*dir) intersects with the region configured as the + // starting point for all vertices. + ///////////////////////////////////////////////////////////////////////// + // Postion in disc + G4double rad = gen_rad_ * std::sqrt(G4UniformRand()); + G4double ang = 2 * G4UniformRand() * pi; + + // Continue assuming that Y is vertical and z drift, + // valid for NEW and NEXT-100 (at least). + // Rotate the disc (origin-point vector) to be perpendicular to dir. + G4ThreeVector point(rad * std::cos(ang), 0., rad * std::sin(ang)); + point.rotate(pi / 2 - dir.angle(point), dir.cross(point)); + + // Now project back to the requested region intersection. + return geom_->ProjectToRegion(region_, point, -dir); +} + + G4double MuonGenerator::GetZenith() const { return fRandomGeneral_->fire()*pi/2; diff --git a/source/generators/MuonGenerator.h b/source/generators/MuonGenerator.h index 36a28572a9..b70482b7c7 100644 --- a/source/generators/MuonGenerator.h +++ b/source/generators/MuonGenerator.h @@ -51,6 +51,8 @@ namespace nexus { void GetDirection(G4ThreeVector& dir, G4double& zenith, G4double& azimuth, G4double& energy, G4double& kinetic_energy, G4double mass); + G4ThreeVector ProjectToVertex(const G4ThreeVector& dir); + G4bool CheckOverlap(const G4ThreeVector& vtx, const G4ThreeVector& dir); /// Load in the Muon Angular/Energy Distribution from CSV file @@ -95,6 +97,8 @@ namespace nexus { std::vector energy_smear_; ///< List of Energy bin smear values G4RandGeneral *fRandomGeneral_; ///< Pointer to the RNG flux distribution + G4double gen_rad_; ///< Radius of disc for generation + }; } // end namespace nexus From 2f088bf4f82baf284a68a71b764604e8e7077dc0 Mon Sep 17 00:00:00 2001 From: Paola Ferrario Date: Tue, 17 Oct 2023 15:34:33 +0200 Subject: [PATCH 22/25] Adapt description to reality of new code --- macros/NEXT100_muons.config.mac | 2 +- macros/NEXT100_muons.init.mac | 6 ++---- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/macros/NEXT100_muons.config.mac b/macros/NEXT100_muons.config.mac index 8fda895e28..48a909462a 100644 --- a/macros/NEXT100_muons.config.mac +++ b/macros/NEXT100_muons.config.mac @@ -1,7 +1,7 @@ ## ---------------------------------------------------------------------------- ## nexus | NEXT100_muons.config.mac ## -## Configuration macro to simulate muons according in the NEXT-100 geometry. +## Configuration macro to simulate muons in the NEXT-100 geometry. ## ## The NEXT Collaboration ## ---------------------------------------------------------------------------- diff --git a/macros/NEXT100_muons.init.mac b/macros/NEXT100_muons.init.mac index 7e6cd645d2..0dd2078534 100644 --- a/macros/NEXT100_muons.init.mac +++ b/macros/NEXT100_muons.init.mac @@ -1,9 +1,7 @@ ## ---------------------------------------------------------------------------- -## nexus | NEXT100_muons_lsc.init.mac +## nexus | NEXT100_muons.init.mac ## -## Initialization macro to simulate muons according to the distribution -## measured at LSC, in the NEXT-100 geometry launching from just outside -## the lead shielding. +## Initialization macro to simulate muons in the NEXT-100 geometry. ## ## The NEXT Collaboration ## ---------------------------------------------------------------------------- From 49a43cb4d5a30a442d8a5eeba915753ef01db39a Mon Sep 17 00:00:00 2001 From: Paola Ferrario Date: Tue, 17 Oct 2023 15:51:15 +0200 Subject: [PATCH 23/25] Fix length of lines for a better readability of the code --- source/geometries/Next100Shielding.cc | 395 ++++++++++++++++---------- 1 file changed, 240 insertions(+), 155 deletions(-) diff --git a/source/geometries/Next100Shielding.cc b/source/geometries/Next100Shielding.cc index 189b04a210..f2d3ff3a95 100644 --- a/source/geometries/Next100Shielding.cc +++ b/source/geometries/Next100Shielding.cc @@ -70,12 +70,15 @@ namespace nexus { { // The shielding is made of two boxes. - // The external one is made of lead and the internal one is made of a mix of stainless steel AISI-316Ti. - // Besides, inside the lead we have the castle steel beams (steel S-275) that support the lead. - // This steel beam structure is composed by tree regions: the lateral beams, the roof, and the - // beam structure above the roof. + // The external one is made of lead and the internal one is made of + // a mix of stainless steel AISI-316Ti. + // Besides, inside the lead we have the castle steel beams (steel S-275) + // that support the lead. + // This steel beam structure is composed by tree regions: + // the lateral beams, the roof, and the beam structure above the roof. - msg_ = new G4GenericMessenger(this, "/Geometry/Next100/", "Control commands of geometry Next100."); + msg_ = new G4GenericMessenger(this, "/Geometry/Next100/", + "Control commands of geometry Next100."); msg_->DeclareProperty("shielding_vis", visibility_, "Shielding Visibility"); msg_->DeclareProperty("shielding_verbosity", verbosity_, "Verbosity"); @@ -91,118 +94,147 @@ namespace nexus { // LEAD BOX /////////// lead_x_ = shield_x_ + 2. * steel_thickness_ + 2. * lead_thickness_; - lead_y_ = shield_y_ + 2. * steel_thickness_ + 2. * lead_thickness_ + beam_thickness_2; + lead_y_ = shield_y_ + 2. * steel_thickness_ + 2. * lead_thickness_ + + beam_thickness_2; lead_z_ = shield_z_ + 2. * steel_thickness_ + 2. * lead_thickness_; - G4Box* lead_box_solid = new G4Box("LEAD_BOX", lead_x_/2., lead_y_/2., lead_z_/2.); + G4Box* lead_box_solid = + new G4Box("LEAD_BOX", lead_x_/2., lead_y_/2., lead_z_/2.); - G4LogicalVolume* lead_box_logic = new G4LogicalVolume(lead_box_solid, - G4NistManager::Instance()->FindOrBuildMaterial("G4_Pb"), - "LEAD_BOX"); + G4LogicalVolume* lead_box_logic = + new G4LogicalVolume(lead_box_solid, + G4NistManager::Instance()->FindOrBuildMaterial("G4_Pb"), + "LEAD_BOX"); this->SetLogicalVolume(lead_box_logic); //STEEL BEAM STRUCTURE // auxiliar positions used in translations - G4double lat_beam_x = shield_x_/2. + steel_thickness_ + lead_thickness_/2.; - G4double front_beam_z = shield_z_/2. + steel_thickness_ + lead_thickness_/2.; - G4double top_beam_y = (shield_y_-beam_thickness_2)/2. + steel_thickness_ + beam_thickness_2 - + lead_thickness_/2.; + G4double lat_beam_x = + shield_x_/2. + steel_thickness_ + lead_thickness_/2.; + G4double front_beam_z = + shield_z_/2. + steel_thickness_ + lead_thickness_/2.; + G4double top_beam_y = + (shield_y_-beam_thickness_2)/2. + steel_thickness_ + beam_thickness_2 + + lead_thickness_/2.; G4double lat_beam_y = -(lead_thickness_ + beam_thickness_2)/2.; - G4double roof_y = (shield_y_-beam_thickness_2)/2. + steel_thickness_ + beam_thickness_2/2.; + G4double roof_y = + (shield_y_-beam_thickness_2)/2. + steel_thickness_ + beam_thickness_2/2.; - G4Box* roof_beam = new G4Box("STRUCT_BEAM", lead_x_/2., beam_thickness_2/2., lead_z_/2.); - G4Box* aux_box = new G4Box("AUX_box", shield_x_/2. + steel_thickness_, beam_thickness_2, - shield_z_/2.+ steel_thickness_); - G4SubtractionSolid* roof_beam_solid = new G4SubtractionSolid("STRUCT_BEAM", roof_beam, aux_box); + G4Box* roof_beam = + new G4Box("STRUCT_BEAM", lead_x_/2., beam_thickness_2/2., lead_z_/2.); + G4Box* aux_box = + new G4Box("AUX_box", shield_x_/2. + steel_thickness_, beam_thickness_2, + shield_z_/2.+ steel_thickness_); + G4SubtractionSolid* roof_beam_solid = + new G4SubtractionSolid("STRUCT_BEAM", roof_beam, aux_box); // vertical bottom beams - G4Box* lat_beam_solid = new G4Box("STRUCT_BEAM", - lead_thickness_/2., - (shield_y_ + 2. * steel_thickness_+ lead_thickness_)/2., - beam_thickness_1/2.); + G4Box* lat_beam_solid = + new G4Box("STRUCT_BEAM", lead_thickness_/2., + (shield_y_ + 2. * steel_thickness_+ lead_thickness_)/2., + beam_thickness_1/2.); - G4Box* front_beam_solid = new G4Box("STRUCT_BEAM", - beam_thickness_2/2., - (shield_y_ + 2. * steel_thickness_+ lead_thickness_)/2., - lead_thickness_/2.); + G4Box* front_beam_solid = + new G4Box("STRUCT_BEAM", beam_thickness_2/2., + (shield_y_ + 2. * steel_thickness_+ lead_thickness_)/2., + lead_thickness_/2.); // horizontal top beams - G4Box* top_xbeam_solid = new G4Box("STRUCT_BEAM", - (shield_x_ + 2.*lead_thickness_ + 2.*steel_thickness_)/2., - lead_thickness_/2., - beam_thickness_1/2.); + G4Box* top_xbeam_solid = + new G4Box("STRUCT_BEAM", + (shield_x_ + 2.*lead_thickness_ + 2.*steel_thickness_)/2., + lead_thickness_/2., beam_thickness_1/2.); - G4Box* top_zbeam_solid = new G4Box("STRUCT_BEAM", - beam_thickness_2/2., - lead_thickness_/2., - (shield_z_ + 2.*lead_thickness_ + 2.*steel_thickness_)/2.); + G4Box* top_zbeam_solid = + new G4Box("STRUCT_BEAM", beam_thickness_2/2., lead_thickness_/2., + (shield_z_ + 2.*lead_thickness_ + 2.*steel_thickness_)/2.); // top beams // x beams G4UnionSolid* struct_solid = - new G4UnionSolid("STEEL_BEAM_STRUCTURE_top1", top_xbeam_solid, top_xbeam_solid, 0, + new G4UnionSolid("STEEL_BEAM_STRUCTURE_top1", + top_xbeam_solid, top_xbeam_solid, 0, G4ThreeVector(0., 0., -roof_z_separation_)); struct_solid = - new G4UnionSolid("STEEL_BEAM_STRUCTURE_top2", struct_solid, top_xbeam_solid, 0, + new G4UnionSolid("STEEL_BEAM_STRUCTURE_top2", + struct_solid, top_xbeam_solid, 0, G4ThreeVector(0., 0., -(roof_z_separation_ + lateral_z_separation_))); struct_solid = - new G4UnionSolid("STEEL_BEAM_STRUCTURE_top3", struct_solid, top_xbeam_solid, 0, + new G4UnionSolid("STEEL_BEAM_STRUCTURE_top3", struct_solid, + top_xbeam_solid, 0, G4ThreeVector(0., 0., -(2*roof_z_separation_ + lateral_z_separation_))); // z beams struct_solid = - new G4UnionSolid("STEEL_BEAM_STRUCTURE_top4", struct_solid, top_zbeam_solid, 0, - G4ThreeVector(-front_x_separation_/2., 0., -(roof_z_separation_+lateral_z_separation_/2.))); + new G4UnionSolid("STEEL_BEAM_STRUCTURE_top4", struct_solid, + top_zbeam_solid, 0, + G4ThreeVector(-front_x_separation_/2., 0., + -(roof_z_separation_+lateral_z_separation_/2.))); struct_solid = - new G4UnionSolid("STEEL_BEAM_STRUCTURE_top5", struct_solid, top_zbeam_solid, 0, - G4ThreeVector(front_x_separation_/2., 0., -(roof_z_separation_+lateral_z_separation_/2.))); + new G4UnionSolid("STEEL_BEAM_STRUCTURE_top5", struct_solid, + top_zbeam_solid, 0, + G4ThreeVector(front_x_separation_/2., 0., + -(roof_z_separation_+lateral_z_separation_/2.))); - G4LogicalVolume* roof_logic = new G4LogicalVolume(roof_beam_solid, - materials::Steel(), "STEEL_BEAM_ROOF"); + G4LogicalVolume* roof_logic = + new G4LogicalVolume(roof_beam_solid, materials::Steel(), + "STEEL_BEAM_ROOF"); - G4LogicalVolume* lat_beam_logic = new G4LogicalVolume(lat_beam_solid, - materials::Steel(), "STEEL_BEAM_STRUCTURE_LAT"); + G4LogicalVolume* lat_beam_logic = + new G4LogicalVolume(lat_beam_solid, materials::Steel(), + "STEEL_BEAM_STRUCTURE_LAT"); - G4LogicalVolume* front_beam_logic = new G4LogicalVolume(front_beam_solid, - materials::Steel(), "STEEL_BEAM_STRUCTURE_FRONT"); + G4LogicalVolume* front_beam_logic = + new G4LogicalVolume(front_beam_solid, materials::Steel(), + "STEEL_BEAM_STRUCTURE_FRONT"); - G4LogicalVolume* struct_logic = new G4LogicalVolume(struct_solid, - materials::Steel(), "STEEL_BEAM_STRUCTURE_TOP"); + G4LogicalVolume* struct_logic = + new G4LogicalVolume(struct_solid, materials::Steel(), + "STEEL_BEAM_STRUCTURE_TOP"); new G4PVPlacement(0, G4ThreeVector(0., top_beam_y, roof_z_separation_ + lateral_z_separation_/2.), - struct_logic, "SHIELDING_STRUCT", lead_box_logic, false, 0, false); + struct_logic, "SHIELDING_STRUCT", lead_box_logic, + false, 0, false); new G4PVPlacement(0, G4ThreeVector(0., roof_y, 0.), - roof_logic, "SHIELDING_STRUCT", lead_box_logic, false, 0, false); + roof_logic, "SHIELDING_STRUCT", lead_box_logic, + false, 0, false); // lateral beams new G4PVPlacement(0, G4ThreeVector(lat_beam_x, lat_beam_y, lateral_z_separation_/2.), - lat_beam_logic, "SHIELDING_STRUCT", lead_box_logic, false, 0, false); + lat_beam_logic, "SHIELDING_STRUCT", lead_box_logic, + false, 0, false); new G4PVPlacement(0, G4ThreeVector(lat_beam_x, lat_beam_y, -lateral_z_separation_/2.), - lat_beam_logic, "SHIELDING_STRUCT", lead_box_logic, false, 0, false); + lat_beam_logic, "SHIELDING_STRUCT", lead_box_logic, + false, 0, false); new G4PVPlacement(0, G4ThreeVector(-lat_beam_x, lat_beam_y, lateral_z_separation_/2.), - lat_beam_logic, "SHIELDING_STRUCT", lead_box_logic, false, 0, false); + lat_beam_logic, "SHIELDING_STRUCT", lead_box_logic, + false, 0, false); new G4PVPlacement(0, G4ThreeVector(-lat_beam_x, lat_beam_y, -lateral_z_separation_/2.), - lat_beam_logic, "SHIELDING_STRUCT", lead_box_logic, false, 0, false); + lat_beam_logic, "SHIELDING_STRUCT", lead_box_logic, + false, 0, false); // front beams new G4PVPlacement(0, G4ThreeVector(-front_x_separation_/2., lat_beam_y, front_beam_z), - front_beam_logic, "SHIELDING_STRUCT", lead_box_logic, false, 0, false); + front_beam_logic, "SHIELDING_STRUCT", lead_box_logic, + false, 0, false); new G4PVPlacement(0, G4ThreeVector(front_x_separation_/2., lat_beam_y, front_beam_z), - front_beam_logic, "SHIELDING_STRUCT", lead_box_logic, false, 0, false); + front_beam_logic, "SHIELDING_STRUCT", lead_box_logic, + false, 0, false); new G4PVPlacement(0, G4ThreeVector(-front_x_separation_/2., lat_beam_y, -front_beam_z), front_beam_logic, "SHIELDING_STRUCT", lead_box_logic, false, 0, false); new G4PVPlacement(0, G4ThreeVector(front_x_separation_/2., lat_beam_y, -front_beam_z), - front_beam_logic,"SHIELDING_STRUCT", lead_box_logic, false, 0, false); + front_beam_logic,"SHIELDING_STRUCT", lead_box_logic, + false, 0, false); // STEEL BOX /////////// @@ -210,7 +242,8 @@ namespace nexus { G4double steel_y = shield_y_ + 2. * steel_thickness_; G4double steel_z = shield_z_ + 2. * steel_thickness_; - G4Box* steel_box_solid = new G4Box("STEEL_BOX", steel_x/2., steel_y/2., steel_z/2.); + G4Box* steel_box_solid = + new G4Box("STEEL_BOX", steel_x/2., steel_y/2., steel_z/2.); G4LogicalVolume* steel_box_logic = new G4LogicalVolume(steel_box_solid, materials::Steel(), "STEEL_BOX"); @@ -220,7 +253,8 @@ namespace nexus { // AIR INSIDE G4Box* air_box_solid = - new G4Box("INNER_AIR", shield_x_/2., shield_y_/2. + steel_thickness_/2., shield_z_/2.); + new G4Box("INNER_AIR", shield_x_/2., shield_y_/2. + steel_thickness_/2., + shield_z_/2.); air_box_logic_ = new G4LogicalVolume(air_box_solid, G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR"), "INNER_AIR"); @@ -243,34 +277,33 @@ namespace nexus { G4double pedestal_support_bottom_thickness = 10. * mm; G4double pedestal_support_top_thickness = 15. * mm; G4double pedestal_support_top_length = 150. * mm; - G4Box* pedestal_support_beam_bottom = new G4Box("PEDESTAL_SUPPORT_BEAM_BOTTOM", - pedestal_x_/2., - lead_thickness_/2., - pedestal_support_bottom_thickness/2.); - - G4Box* pedestal_support_beam_top = new G4Box("PEDESTAL_SUPPORT_BEAM_TOP", - pedestal_top_x_/2., - pedestal_support_top_thickness/2., - pedestal_support_top_length/2.); - - G4Box* pedestal_beam_front = new G4Box("PEDESTAL_BEAM_FRONT", - pedestal_x_/2., - lead_thickness_/2., - pedestal_front_beam_thickness_/2.); - - G4Box* pedestal_beam_lateral = new G4Box("PEDESTAL_BEAM_LATERAL", - pedestal_lateral_beam_thickness_/2., - lead_thickness_/2., - pedestal_lateral_length_/2.); + G4Box* pedestal_support_beam_bottom = + new G4Box("PEDESTAL_SUPPORT_BEAM_BOTTOM", + pedestal_x_/2., lead_thickness_/2., + pedestal_support_bottom_thickness/2.); + + G4Box* pedestal_support_beam_top = + new G4Box("PEDESTAL_SUPPORT_BEAM_TOP", pedestal_top_x_/2., + pedestal_support_top_thickness/2., + pedestal_support_top_length/2.); + + G4Box* pedestal_beam_front = + new G4Box("PEDESTAL_BEAM_FRONT", pedestal_x_/2., lead_thickness_/2., + pedestal_front_beam_thickness_/2.); + + G4Box* pedestal_beam_lateral = + new G4Box("PEDESTAL_BEAM_LATERAL", pedestal_lateral_beam_thickness_/2., + lead_thickness_/2., pedestal_lateral_length_/2.); G4Box* pedestal_roof_ = new G4Box("PEDESTAL_ROOF", pedestal_x_/2. + pedestal_lateral_beam_thickness_, - pedestal_lateral_beam_thickness_/2., pedestal_lateral_length_/2.); + pedestal_lateral_beam_thickness_/2., + pedestal_lateral_length_/2.); G4Box* ped_aux_box = new G4Box("AUX_box", pedestal_x_/2. + pedestal_lateral_beam_thickness_ - pedestal_roof_thickness_, - pedestal_lateral_beam_thickness_, - pedestal_lateral_length_/2.-pedestal_roof_thickness_); + pedestal_lateral_beam_thickness_, + pedestal_lateral_length_/2.-pedestal_roof_thickness_); G4SubtractionSolid* pedestal_roof = new G4SubtractionSolid("PEDESTAL_ROOF", pedestal_roof_, ped_aux_box); @@ -297,10 +330,10 @@ namespace nexus { G4double pedestal_x_pos = pedestal_x_/2. + pedestal_lateral_beam_thickness_/2.; G4double pedestal_y_pos = -lead_y_/2. + lead_thickness_/2.; - G4double pedestal_top_y_pos = -(shield_y_/2. + steel_thickness_/2.) - + pedestal_support_top_thickness/2.; //caution: it is refered to air-box - G4double pedestal_roof_y_pos = -(shield_y_/2. + steel_thickness_/2.) - + pedestal_lateral_beam_thickness_/2.; //caution: it is refered to air-box + G4double pedestal_top_y_pos = -(shield_y_/2. + steel_thickness_/2.) + + pedestal_support_top_thickness/2.; //caution: it is refered to air-box + G4double pedestal_roof_y_pos = -(shield_y_/2. + steel_thickness_/2.) + + pedestal_lateral_beam_thickness_/2.; //caution: it is refered to air-box new G4PVPlacement(0, G4ThreeVector(0., pedestal_y_pos, support_beam_dist_/2.), @@ -378,8 +411,9 @@ namespace nexus { shield_diag / 2. + ext_offset, 1. * mm, G4ThreeVector(0.,0.,0.), 0); - steel_gen_ = new BoxPointSampler(shield_x_, shield_y_, shield_z_, steel_thickness_, - G4ThreeVector(0., -beam_thickness_2/2., 0.), 0); + steel_gen_ = + new BoxPointSampler(shield_x_, shield_y_, shield_z_, steel_thickness_, + G4ThreeVector(0., -beam_thickness_2/2., 0.), 0); inner_air_gen_ = new BoxPointSampler(shield_x_, shield_y_, shield_z_, 0, G4ThreeVector(0., 0., 0.), 0); @@ -387,7 +421,8 @@ namespace nexus { // STEEL STRUCTURE GENERATORS lat_roof_gen_ = - new BoxPointSampler(lead_thickness_, beam_thickness_2, shield_z_ + 2.*steel_thickness_, 0., + new BoxPointSampler(lead_thickness_, beam_thickness_2, + shield_z_ + 2.*steel_thickness_, 0., G4ThreeVector(0., roof_y, 0.), 0); front_roof_gen_ = @@ -395,24 +430,29 @@ namespace nexus { G4ThreeVector(0., roof_y, 0.), 0); struct_x_gen_ = - new BoxPointSampler(shield_x_ + 2.*lead_thickness_ + 2.*steel_thickness_, lead_thickness_, - beam_thickness_1, 0., + new BoxPointSampler(shield_x_ + 2.*lead_thickness_ + 2.*steel_thickness_, + lead_thickness_, beam_thickness_1, 0., G4ThreeVector(0., top_beam_y, roof_z_separation_+lateral_z_separation_/2.), 0); struct_z_gen_ = new BoxPointSampler(beam_thickness_2, lead_thickness_, - shield_z_ + 2.*lead_thickness_ + 2.*steel_thickness_, 0., + shield_z_ + 2.*lead_thickness_ + 2.*steel_thickness_, + 0., G4ThreeVector(-front_x_separation_/2., top_beam_y, 0.), 0); lat_beam_gen_ = - new BoxPointSampler(lead_thickness_, shield_y_ + 2. * steel_thickness_+ lead_thickness_, + new BoxPointSampler(lead_thickness_, + shield_y_ + 2. * steel_thickness_+ lead_thickness_, beam_thickness_1, 0., - G4ThreeVector(lat_beam_x, -lead_thickness_/2., lateral_z_separation_/2.), 0); + G4ThreeVector(lat_beam_x, -lead_thickness_/2., + lateral_z_separation_/2.), 0); front_beam_gen_ = - new BoxPointSampler(beam_thickness_2, shield_y_ + 2. * steel_thickness_+ lead_thickness_, + new BoxPointSampler(beam_thickness_2, + shield_y_ + 2. * steel_thickness_+ lead_thickness_, lead_thickness_, 0., - G4ThreeVector(-front_x_separation_/2., -lead_thickness_/2., front_beam_z), 0); + G4ThreeVector(-front_x_separation_/2., + -lead_thickness_/2., front_beam_z), 0); // Compute relative volumes G4double roof_vol = roof_beam_solid->GetCubicVolume(); @@ -424,35 +464,43 @@ namespace nexus { perc_front_roof_vol_ = 2*(lead_x_*beam_thickness_2*lead_thickness_)/roof_vol; perc_top_struct_vol_ = struct_top_vol /total_vol; - G4double struc_beam_x_vol = (shield_x_ + 2.*lead_thickness_ - + 2.*steel_thickness_)*lead_thickness_*beam_thickness_1; - perc_struc_x_vol_ = 4*struc_beam_x_vol/struct_top_vol; + G4double struc_beam_x_vol = + (shield_x_ + 2.*lead_thickness_ + + 2.*steel_thickness_)*lead_thickness_*beam_thickness_1; + perc_struc_x_vol_ = 4*struc_beam_x_vol/struct_top_vol; // PEDESTAL GENERATORS ped_support_bottom_gen_ = - new BoxPointSampler(pedestal_x_, lead_thickness_, pedestal_support_bottom_thickness, 0., + new BoxPointSampler(pedestal_x_, lead_thickness_, + pedestal_support_bottom_thickness, 0., G4ThreeVector(0., pedestal_y_pos, support_beam_dist_/2.), 0); G4double ped_gen_y_ = -lead_y_/2. + lead_thickness_ + pedestal_support_top_thickness/2.; ped_support_top_gen_ = - new BoxPointSampler(pedestal_top_x_, pedestal_support_top_thickness, pedestal_support_top_length, 0., + new BoxPointSampler(pedestal_top_x_, pedestal_support_top_thickness, + pedestal_support_top_length, 0., G4ThreeVector(0., ped_gen_y_, support_beam_dist_/2.), 0); ped_front_gen_ = - new BoxPointSampler(pedestal_x_, lead_thickness_, pedestal_front_beam_thickness_, 0., + new BoxPointSampler(pedestal_x_, lead_thickness_, + pedestal_front_beam_thickness_, 0., G4ThreeVector(0., pedestal_y_pos, support_beam_dist_/2. + support_front_dist_), 0); ped_lateral_gen_ = - new BoxPointSampler(pedestal_lateral_beam_thickness_, lead_thickness_, pedestal_lateral_length_, 0., + new BoxPointSampler(pedestal_lateral_beam_thickness_, lead_thickness_, + pedestal_lateral_length_, 0., G4ThreeVector(pedestal_x_pos, pedestal_y_pos, 0.), 0); // pedestal roof G4double ped_roof_gen_x = pedestal_top_x_/2. + pedestal_roof_thickness_/2.; - G4double ped_roof_gen_y = -lead_y_/2. + lead_thickness_ + pedestal_lateral_beam_thickness_/2.; - G4double ped_roof_gen_z = pedestal_lateral_length_/2. - pedestal_roof_thickness_/2.; + G4double ped_roof_gen_y = -lead_y_/2. + lead_thickness_ + + pedestal_lateral_beam_thickness_/2.; + G4double ped_roof_gen_z = pedestal_lateral_length_/2. - + pedestal_roof_thickness_/2.; ped_roof_lat_gen_ = - new BoxPointSampler(pedestal_roof_thickness_, pedestal_lateral_beam_thickness_, + new BoxPointSampler(pedestal_roof_thickness_, + pedestal_lateral_beam_thickness_, pedestal_lateral_length_, 0., G4ThreeVector(ped_roof_gen_x, ped_roof_gen_y, 0.), 0); @@ -467,8 +515,10 @@ namespace nexus { G4double ped_lateral_vol = pedestal_beam_lateral ->GetCubicVolume(); G4double ped_roof_vol = pedestal_roof ->GetCubicVolume(); - G4double ped_total_vol = ped_roof_vol + 2*(ped_support_bottom_vol + ped_support_top_vol - + ped_front_vol + ped_lateral_vol); + G4double ped_total_vol = ped_roof_vol + 2*(ped_support_bottom_vol + + ped_support_top_vol + + ped_front_vol + + ped_lateral_vol); perc_ped_bottom_vol_ = 2*ped_support_bottom_vol/ped_total_vol; perc_ped_top_vol_ = 2*ped_support_top_vol /ped_total_vol; @@ -479,36 +529,47 @@ namespace nexus { // BUBBLE SEAL bubble_seal_front_gen_ = - new BoxPointSampler(pedestal_x_ + 2*pedestal_lateral_beam_thickness_, bubble_seal_thickness_, + new BoxPointSampler(pedestal_x_ + 2*pedestal_lateral_beam_thickness_, + bubble_seal_thickness_, bubble_seal_thickness_, 0., G4ThreeVector(0., ped_gen_y_, 0.), 0); bubble_seal_lateral_gen_ = - new BoxPointSampler(bubble_seal_thickness_, bubble_seal_thickness_, pedestal_lateral_length_, 0., + new BoxPointSampler(bubble_seal_thickness_, bubble_seal_thickness_, + pedestal_lateral_length_, 0., G4ThreeVector(0., ped_gen_y_, 0.), 0); - G4double bubble_front_vol = 2*(pedestal_x_ + 2*pedestal_lateral_beam_thickness_) - *(bubble_seal_thickness_)*(bubble_seal_thickness_); - G4double bubble_lateral_vol = 2*(bubble_seal_thickness_) - *(bubble_seal_thickness_)*(pedestal_lateral_length_); + G4double bubble_front_vol = + 2*(pedestal_x_ + 2*pedestal_lateral_beam_thickness_) + *(bubble_seal_thickness_)*(bubble_seal_thickness_); + G4double bubble_lateral_vol = + 2*(bubble_seal_thickness_) + *(bubble_seal_thickness_)*(pedestal_lateral_length_); - perc_bubble_front_vol_ = bubble_front_vol /(bubble_front_vol + bubble_lateral_vol); - perc_buble_lateral_vol_ = bubble_lateral_vol/(bubble_front_vol + bubble_lateral_vol); + perc_bubble_front_vol_ = + bubble_front_vol /(bubble_front_vol + bubble_lateral_vol); + perc_buble_lateral_vol_= + bubble_lateral_vol/(bubble_front_vol + bubble_lateral_vol); // EDPM SEAL edpm_seal_front_gen_ = - new BoxPointSampler(edpm_seal_thickness_, shield_y_, edpm_seal_thickness_, 0., + new BoxPointSampler(edpm_seal_thickness_, shield_y_, + edpm_seal_thickness_, 0., G4ThreeVector(0., -(beam_thickness_2+steel_thickness_)/2., 0.), 0); edpm_seal_lateral_gen_ = - new BoxPointSampler(edpm_seal_thickness_, edpm_seal_thickness_, shield_z_, 0., + new BoxPointSampler(edpm_seal_thickness_, edpm_seal_thickness_, + shield_z_, 0., G4ThreeVector(0., -(beam_thickness_2+steel_thickness_)/2., 0.), 0); - G4double edpm_front_vol = 2*edpm_seal_thickness_ * shield_y_ * edpm_seal_thickness_; - G4double edpm_lateral_vol = edpm_seal_thickness_ * edpm_seal_thickness_ * shield_z_; + G4double edpm_front_vol = + 2*edpm_seal_thickness_ * shield_y_ * edpm_seal_thickness_; + G4double edpm_lateral_vol = + edpm_seal_thickness_ * edpm_seal_thickness_ * shield_z_; - perc_edpm_front_vol_ = edpm_front_vol /(edpm_front_vol + edpm_lateral_vol); - perc_edpm_lateral_vol_ = edpm_lateral_vol/(edpm_front_vol + edpm_lateral_vol); + perc_edpm_front_vol_ = edpm_front_vol /(edpm_front_vol + edpm_lateral_vol); + perc_edpm_lateral_vol_ = + edpm_lateral_vol/(edpm_front_vol + edpm_lateral_vol); if (verbosity_){ @@ -583,20 +644,22 @@ namespace nexus { if (region == "SHIELDING_LEAD") { G4VPhysicalVolume *VertexVolume; do { - vertex = lead_gen_->GenerateVertex("WHOLE_VOL"); - G4ThreeVector glob_vtx(vertex); - glob_vtx = glob_vtx + G4ThreeVector(0, 0, -GetELzCoord()); - VertexVolume = geom_navigator_->LocateGlobalPointAndSetup(glob_vtx, 0, false); + vertex = lead_gen_->GenerateVertex("WHOLE_VOL"); + G4ThreeVector glob_vtx(vertex); + glob_vtx = glob_vtx + G4ThreeVector(0, 0, -GetELzCoord()); + VertexVolume = geom_navigator_->LocateGlobalPointAndSetup(glob_vtx, + 0, false); } while (VertexVolume->GetName() != "LEAD_BOX"); } else if (region == "SHIELDING_STEEL") { G4VPhysicalVolume *VertexVolume; do { - vertex = steel_gen_->GenerateVertex("WHOLE_VOL"); - G4ThreeVector glob_vtx(vertex); - glob_vtx = glob_vtx + G4ThreeVector(0, 0, -GetELzCoord()); - VertexVolume = geom_navigator_->LocateGlobalPointAndSetup(glob_vtx, 0, false); + vertex = steel_gen_->GenerateVertex("WHOLE_VOL"); + G4ThreeVector glob_vtx(vertex); + glob_vtx = glob_vtx + G4ThreeVector(0, 0, -GetELzCoord()); + VertexVolume = geom_navigator_->LocateGlobalPointAndSetup(glob_vtx, + 0, false); } while (VertexVolume->GetName() != "STEEL_BOX"); } @@ -606,7 +669,8 @@ namespace nexus { vertex = inner_air_gen_->GenerateVertex("INSIDE"); G4ThreeVector glob_vtx(vertex); glob_vtx = glob_vtx + G4ThreeVector(0, 0, -GetELzCoord()); - VertexVolume = geom_navigator_->LocateGlobalPointAndSetup(glob_vtx, 0, false); + VertexVolume = geom_navigator_->LocateGlobalPointAndSetup(glob_vtx, + 0, false); } while (VertexVolume->GetName() != "INNER_AIR"); } @@ -621,19 +685,25 @@ namespace nexus { if (G4UniformRand() < perc_front_roof_vol_){ vertex = front_roof_gen_->GenerateVertex("INSIDE"); if (G4UniformRand() < 0.5) { - vertex.setZ(vertex.z() + (shield_z_/2.+steel_thickness_+lead_thickness_/2.)); + vertex.setZ(vertex.z() + + (shield_z_/2.+steel_thickness_+lead_thickness_/2.)); } else{ - vertex.setZ(vertex.z() - (shield_z_/2.+steel_thickness_+lead_thickness_/2.)); + vertex.setZ(vertex.z() - + (shield_z_/2.+steel_thickness_+lead_thickness_/2.)); } } else{ vertex = lat_roof_gen_->GenerateVertex("INSIDE"); if (G4UniformRand() < 0.5) { - vertex.setX(vertex.x() + (shield_x_/2.+ steel_thickness_ + lead_thickness_/2.)); + vertex.setX(vertex.x() + + (shield_x_/2.+ steel_thickness_ + + lead_thickness_/2.)); } else{ - vertex.setX(vertex.x() - (shield_x_/2.+ steel_thickness_ + lead_thickness_/2.)); + vertex.setX(vertex.x() - + (shield_x_/2.+ steel_thickness_ + + lead_thickness_/2.)); } } } @@ -647,10 +717,12 @@ namespace nexus { vertex.setZ(vertex.z()-roof_z_separation_); } else if (rand_beam == 2) { - vertex.setZ(vertex.z()-(roof_z_separation_+lateral_z_separation_)); + vertex.setZ(vertex.z()-(roof_z_separation_ + + lateral_z_separation_)); } else if (rand_beam == 3) { - vertex.setZ(vertex.z()-(2*roof_z_separation_+lateral_z_separation_)); + vertex.setZ(vertex.z()-(2*roof_z_separation_ + + lateral_z_separation_)); } } else { @@ -670,10 +742,12 @@ namespace nexus { vertex.setZ(vertex.z() - lateral_z_separation_); } else if (rand_beam ==2){ - vertex.setX(vertex.x() - (shield_x_ + 2*steel_thickness_ + lead_thickness_)); + vertex.setX(vertex.x() - (shield_x_ + 2*steel_thickness_ + + lead_thickness_)); } else if (rand_beam ==3){ - vertex.setX(vertex.x() - (shield_x_ + 2*steel_thickness_ + lead_thickness_)); + vertex.setX(vertex.x() - (shield_x_ + 2*steel_thickness_ + + lead_thickness_)); vertex.setZ(vertex.z() - lateral_z_separation_); } } @@ -684,11 +758,13 @@ namespace nexus { vertex.setX(vertex.x() + front_x_separation_); } else if (rand_beam ==2){ - vertex.setZ(vertex.z() - (shield_z_+2*steel_thickness_+lead_thickness_)); + vertex.setZ(vertex.z() - (shield_z_+2*steel_thickness_ + + lead_thickness_)); } else if (rand_beam ==3){ vertex.setX(vertex.x() + front_x_separation_); - vertex.setZ(vertex.z() - (shield_z_+2*steel_thickness_+lead_thickness_)); + vertex.setZ(vertex.z() - (shield_z_+2*steel_thickness_ + + lead_thickness_)); } } } @@ -709,30 +785,35 @@ namespace nexus { vertex.setZ(vertex.z() - support_beam_dist_); } } - else if (rand < (perc_ped_bottom_vol_ + perc_ped_top_vol_ + perc_ped_front_vol_)){ //FRONT BEAM + else if (rand < (perc_ped_bottom_vol_ + perc_ped_top_vol_ + + perc_ped_front_vol_)){ //FRONT BEAM vertex = ped_front_gen_->GenerateVertex("INSIDE"); if (G4UniformRand() < 0.5) { - vertex.setZ(vertex.z() - (2.*support_front_dist_ + support_beam_dist_)); + vertex.setZ(vertex.z() - (2.*support_front_dist_ + + support_beam_dist_)); } } - else if (rand < (perc_ped_bottom_vol_ + perc_ped_top_vol_ + perc_ped_front_vol_ - + perc_ped_lateral_vol_)){ // LATERAL BEAM + else if (rand < (perc_ped_bottom_vol_ + perc_ped_top_vol_ + + perc_ped_front_vol_ + perc_ped_lateral_vol_)){ // LATERAL BEAM vertex = ped_lateral_gen_->GenerateVertex("INSIDE"); if (G4UniformRand() < 0.5) { - vertex.setX(vertex.x() - (pedestal_x_ + pedestal_lateral_beam_thickness_)); + vertex.setX(vertex.x() - (pedestal_x_ + + pedestal_lateral_beam_thickness_)); } } else { // ROOF if (G4UniformRand() < 0.5) { vertex = ped_roof_lat_gen_->GenerateVertex("INSIDE"); if (G4UniformRand() < 0.5){ - vertex.setX(vertex.x() - pedestal_top_x_ - pedestal_roof_thickness_); + vertex.setX(vertex.x() - pedestal_top_x_ - + pedestal_roof_thickness_); } } else{ vertex = ped_roof_front_gen_->GenerateVertex("INSIDE"); if (G4UniformRand() < 0.5){ - vertex.setZ(vertex.z() - pedestal_lateral_length_ + pedestal_roof_thickness_); + vertex.setZ(vertex.z() - pedestal_lateral_length_ + + pedestal_roof_thickness_); } } } @@ -746,22 +827,26 @@ namespace nexus { vertex = bubble_seal_front_gen_->GenerateVertex("INSIDE"); if (G4UniformRand() < 0.5){ vertex.setZ(vertex.z() + (support_beam_dist_/2. + support_front_dist_ - + pedestal_front_beam_thickness_/2. + bubble_seal_thickness_/2.)); + + pedestal_front_beam_thickness_/2. + + bubble_seal_thickness_/2.)); } else{ vertex.setZ(vertex.z() - (support_beam_dist_/2. + support_front_dist_ - + pedestal_front_beam_thickness_/2. + bubble_seal_thickness_/2.)); + + pedestal_front_beam_thickness_/2. + + bubble_seal_thickness_/2.)); } } else{ // lateral vertex = bubble_seal_lateral_gen_->GenerateVertex("INSIDE"); if (G4UniformRand() < 0.5){ - vertex.setX(vertex.x() + (pedestal_x_/2. + pedestal_lateral_beam_thickness_ - + bubble_seal_thickness_/2.)); + vertex.setX(vertex.x() + (pedestal_x_/2. + + pedestal_lateral_beam_thickness_ + + bubble_seal_thickness_/2.)); } else{ - vertex.setX(vertex.x() - (pedestal_x_/2. + pedestal_lateral_beam_thickness_ - + bubble_seal_thickness_/2.)); + vertex.setX(vertex.x() - (pedestal_x_/2. + + pedestal_lateral_beam_thickness_ + + bubble_seal_thickness_/2.)); } } } From e4e1622ac7a02ee51341b3e784a3b102adf05d48 Mon Sep 17 00:00:00 2001 From: Paola Ferrario Date: Tue, 17 Oct 2023 15:51:57 +0200 Subject: [PATCH 24/25] Change variable name to avoid confusion con CLHEP::rad --- source/generators/MuonGenerator.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/generators/MuonGenerator.cc b/source/generators/MuonGenerator.cc index 158f474e1d..8f2aaabebe 100644 --- a/source/generators/MuonGenerator.cc +++ b/source/generators/MuonGenerator.cc @@ -353,13 +353,13 @@ G4ThreeVector MuonGenerator::ProjectToVertex(const G4ThreeVector& dir) // starting point for all vertices. ///////////////////////////////////////////////////////////////////////// // Postion in disc - G4double rad = gen_rad_ * std::sqrt(G4UniformRand()); + G4double radius = gen_rad_ * std::sqrt(G4UniformRand()); G4double ang = 2 * G4UniformRand() * pi; // Continue assuming that Y is vertical and z drift, // valid for NEW and NEXT-100 (at least). // Rotate the disc (origin-point vector) to be perpendicular to dir. - G4ThreeVector point(rad * std::cos(ang), 0., rad * std::sin(ang)); + G4ThreeVector point(radius * std::cos(ang), 0., radius * std::sin(ang)); point.rotate(pi / 2 - dir.angle(point), dir.cross(point)); // Now project back to the requested region intersection. From 224461437c2666372f2857762ee7bb29f3bf564b Mon Sep 17 00:00:00 2001 From: Paola Ferrario Date: Tue, 17 Oct 2023 16:08:47 +0200 Subject: [PATCH 25/25] Revert the change in EXTERNAL gen because we don't see the point --- source/geometries/Next100Shielding.cc | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/source/geometries/Next100Shielding.cc b/source/geometries/Next100Shielding.cc index f2d3ff3a95..0d5b05e815 100644 --- a/source/geometries/Next100Shielding.cc +++ b/source/geometries/Next100Shielding.cc @@ -403,12 +403,10 @@ namespace nexus { // Only shooting from the innest 5 cm. lead_gen_ = new BoxPointSampler(steel_x, steel_y, steel_z, 5.*cm, G4ThreeVector(0., 0., 0.), 0); - - G4double shield_diag = std::sqrt(lead_x_*lead_x_ + lead_y_*lead_y_ + lead_z_*lead_z_); G4double ext_offset = 1. * cm; - external_gen_ = new BoxPointSampler(shield_diag / 2. + ext_offset, - shield_diag / 2. + ext_offset, - shield_diag / 2. + ext_offset, + external_gen_ = new BoxPointSampler(lead_x_ + ext_offset, + lead_y_ + ext_offset, + lead_z_ + ext_offset, 1. * mm, G4ThreeVector(0.,0.,0.), 0); steel_gen_ =